A combination of convergent extension and differential adhesion explains the shapes of elongating gastruloids¶

This notebook was created in Jupyter lab. The required Python packages are stored in the environment file environment_gastrushape.yml. To import this environment, see the instructions in step 0.1 and 0.2, but replace the name of the yml-file with environment_gastrushape.yml and the environment name with gastrushape.

Import the required libraries¶

In [2]:
# The basic libraries
import numpy as np 
import pandas as pd
import os
import math
import scipy

# To plot data
import matplotlib
from matplotlib import pyplot as plt
import matplotlib.image as mpimg
matplotlib.use('Agg')
import sys

# To read in images
import imageio
# To process images
import cv2

# To read ROIs generated in ImageJ
from read_roi import read_roi_file
from read_roi import read_roi_zip

# To fit ellipse
from skimage.measure import EllipseModel
from matplotlib.patches import Ellipse

LOCO-EFA of gastruloid shapes¶

Calculate the LOCO-EFA values of gastruloids shapes¶

We obtained the LOCO-EFA values from gastruloid shapes that were extracted from Nikon ND2 files. Our ND2 files often contain images from multiple gastruloids in one field of view, of which we recorded multiple fluorescent channels, multiple positions in Z (z-stack), and multiple positions in XY. We developed a segmentation tool that deals with this format and extracts outlines from individual gastruloid shapes.

To run this tool, follow below steps. The tool may also be used for image files with a different extension, such as .jpg, .png, .tif or .tiff. If the image file contains a bright-field image, segmentation will likely not be performed correctly. The signal in the image needs to be limited to the gastruloid and should not be visible in the background. Therefore, the image file should only contain fluorescent channels, ideally one that is homogeneously showing the gastruloids shape, for example by using DAPI to visualize all nuclei.

NB: Make sure that your images contain unique gastruloids, or that you remove 'doubles' from the analysis later on. In case you imaged multiple gastruloids and took several images with overlapping areas, there's a chance that two different images include the same gastruloid. Make sure you remove 'doubles' from the LOCO-EFA analysis to not have them in the final plot. We carefully checked for 'doubles' and did not include them in the plot.

0. Create the environment from the environment.yml file¶

All the required files are stored in the Interface folder. Navigate in your terminal to this folder using the cd tool. In this folder you will find the following files:

  • environment.yml
  • ciha.py
  • gui-segmentation.py

How to run the segmentation tool gui-segmentation.py?

  1. You first have to create the environment from the environment.yml file by running
      conda env create -f environment.yml
    The name of the environment is ciha. You can check if it's in the conda environments list by running
      conda env list

  2. You then activate the environment by running
      source activate ciha

  3. You run the segmentation tool by running
      python gui-segmentation.py
    The segmentation tool imports ciha.py, which defines several functions.

The next time you want to use the segmentation tool, you only need to navigate to the interface folder and run step 2 and 3.

1. Use our segmentation tool to create outlines from fluorescently labelled gastruloids¶

Explanation of buttons (in order of use)
File Open: opens a file dialog window to let you select a folder. If the folder contains multiple image files, all will be recognized when they have the extension .nd2, .png, .jpg, .tif, or .tiff.

Run segmentation: runs the segmentation on the first image in the folder with the standard parameters (threshold value = none, minimal distance = 70, minimal size = 100, Gaussian smoothing = 5, opening for maxima = 5 [pixels?]). These parameters can be adjusted if desired.

For each image position (XY), save the outlines:

File Save Outlines: opens a file dialog window. Select the folder in which the outlines should be stored and provide a name for the file including the extension (e.g., .txt or .csv).

Save Outlines: saves the outlines that were found after segmentation in the selected folder with the set name.

Next xy: If an image file contains multiple XY positions, you may click Next xy to go to the next position and repeat the segmentation procedure, the selection of the folder and naming the file, and saving the outlines.

Next File: if the folder contains multiple images, you select the next image by clicking this button. The segmentation tool keeps track of the image number and the position number in the boxes Current Image and Current xy, respectively.

Run Fourier Analysis: runs a Fourier analysis on the outlines. The calculated Fourier components can be stored using the buttons File Save Fourier to select the path and file name, and Save Fourier Components to store the components in the selected path.

The segmentation tool stores all the outlines of the detected gastruloids in a single .csv or .txt file. To generate from outlines LOCO-EFA values that describe the shape, each outline should be stored as a separate .csv file. To do so, you may go over the following part of this jupyter notebook that stores the outlines separately.

2. Store each outline into a separate csv file¶

In [ ]:
# Import the required libraries
from numpy import array, savetxt
In [ ]:
# Select the txt file (a file dialog will open)
try:
    from tkinter import Tk
    from tkFileDialog import askopenfilenames
except:
    from tkinter import Tk
    from tkinter import filedialog

Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
filenames = filedialog.askopenfilenames() # show an "Open" dialog box and return the path to the selected file
filename = ''.join(filenames)
print(filename)
In [ ]:
# Read in file, and get the number of outlines listed in the file
reader = open(filename,'r')
try:
    #print(reader.read())
    value_lines = list(reader)
finally:
    reader.close()
size_val = len(value_lines)

# Get the number of outlines listed in the file and check, if not correct, adjust it
vals = value_lines[size_val-1].strip("\n")
vals2 = vals.split("\t")
n_outlines = int(vals2[0][-1])+1
print('Number of outlines in this file are:', n_outlines)
#n_outlines = 3
#print(n_outlines)
In [ ]:
# Create a new folder to store the separate outlines in
# Get the file path (file_path) and file name (report_name) without the extension
filename_split = filename.split("/")
filename_split_length = len(filename_split)
file_path = (filename[0:(len(filename)-len(filename_split[filename_split_length-1]))])
print(file_path)

report_nametxt = filename_split[filename_split_length-1]
report_name = report_nametxt[0:(len(report_nametxt)-4)]
print(report_name)

# Create a new folder to store the header and the separate outline files in
newpath = '%sFolder_%s' % (file_path,report_name)
print(newpath)
if not os.path.exists(newpath):
    os.makedirs(newpath)

In case you have many imaging files, you will also have many gastruloid outlines. To save outlines with a file number starting from a certain number (to eventually have outline_1 to outline_100 instead of 10 outline_1 files), change count_start to the desired number. Keep a record on this starting number for each new XY, and subsequently each new imaging file. Note that the code in the cell below doesn't work yet for txt or csv files that contain more than 10 outlines. To handle those, first split the files in two (0-9 and 10-19) or more (20-29, etc.).
>> currently the code only reads out the last value of the csv or txt file (numbers 0-9) to determine the number of outlines that are listed in this file.

In [ ]:
# To save outlines with a file number, starting from a certain number, change count_start to the desired number.
count_start = 10;
In [ ]:
# For each outline (j), list the x and y values in separate arrays and write the arrays to a csv file
for j in range(0,n_outlines):
    x_array = np.array([])
    y_array = np.array([])
    
    for i in range(1,size_val):
       # print('i= %i and j= %j',(i,j))
        value_lines_stripped = value_lines[i].strip("\n")
        value_lines_split = value_lines_stripped.split("\t")
    
        if j == (n_outlines-1) and i == (size_val-1): # j is the last outline and i is the last written line 
            # >> append array, save array, save header
            xval = float(value_lines_split[1])
            xvalr = round(xval)
            yval = float(value_lines_split[2])
            yvalr = round(yval)
            x_array = np.append(x_array,xvalr)
            y_array = np.append(y_array,yvalr)
            header_outline = value_lines[0].strip("\n")
            np.savetxt('%s/Outline_%i.csv' % (newpath, j+count_start), np.array([x_array, y_array]).T, fmt='%i', delimiter=',')#, header = header_cont)
            np.savetxt('%s/Header_%s.csv' % (newpath,report_name), [0], header = header_outline)
        elif int(value_lines_split[0][-1]) == j: # if j is not yet 9, and line i belongs to j >> append array
            xval = float(value_lines_split[1])
            xvalr = round(xval)
            yval = float(value_lines_split[2])
            yvalr = round(yval)
            x_array = np.append(x_array,xvalr)
            y_array = np.append(y_array,yvalr)
        elif int(value_lines_split[0][-1]) > j: # if j is not yet 9, and line i is larger than j >> save text
            np.savetxt('%s/Outline_%i.csv' % (newpath, j+count_start), np.array([x_array, y_array]).T, fmt='%i', delimiter=',')#, header = header_cont  
            i = size_val

It can be helpful to visualize the gastruloid shape from the outline file. The following code allows for that.

In [4]:
outline_data = pd.read_csv('Outlines/outlines138', header=None)
n_pixels_outline = round(outline_data.size / 2)
outline_matrix = np.zeros((1024,1024))
for i in range(n_pixels_outline):
    read_x = round(outline_data.loc[i,0])
    read_y = round(outline_data.loc[i,1])
    outline_matrix[read_x,read_y] = 1

# Fill the outline image
filled_binary = scipy.ndimage.binary_fill_holes(outline_matrix)

# Plot shape
%matplotlib inline
plt.imshow(filled_binary, cmap = 'gray')
Out[4]:
<matplotlib.image.AxesImage at 0x7f9ea16c3f10>

3. Calculate for each outline file the LOCO-EFA values¶

Collect all the outline files into one folder.

The LOCO-EFA method was created by Yara E Sánchez-Corrales, Matthew Hartley, Jop van Rooij, Athanasius F M Marée, Verônica A Grieneisen (Development. 2018 Mar 20;145(6):dev156778. doi: 10.1242/dev.156778.) LOCO-EFA stands for 'Lobe contribution elliptical Fourier analysis (LOCO-EFA)'. This method generates from digitalised two-dimensional cell outlines meaningful descriptors that can be directly matched to morphological features.

Link to paper: https://pubmed.ncbi.nlm.nih.gov/29444894/
Their (example) code can be found here: https://bitbucket.org/mareelab/loco_efa/src/master/src/loco-efa_example.c

We slightly adjusted the code so that the L-values found for each outline were stored in a text file in a path that we set in line 17 in our .c file. When using our file, you should adjust this path to a path on your pc. Name of our file: loco-efa_example_GS.c
Below's code worked on our MacBook Pro from 2017.
First, the code needs to be compiled (C is a mid-level language and it needs a compiler to convert it into an executable code so that the program can be run on our machine).

  1. Navigate to the folder that includes the file 'loco-efa_example_GS.c', using cd and pwd
  2. Run the following command:
      cc -Wall -o loco-efa_example_GS loco-efa_example_GS.c -lm
  3. If you get the error 'xcrun: error: invalid active developer path (/Library/Developer/CommandLineTools), missing xcrun at: /Library/Developer/CommandLineTools/usr/bin/xcrun' you need to install the command line tools for Mac.
    <br>   a. Check the macOS version of your system and check online the release date of this version. 
    <br>   b. Go to https://developer.apple.com/download/all/?q=command%20line%20tools and download the command line tools for Xcode that were released before the release date of your latest macOS version.
     c. Install the command line tools
  4. Repeat step 1 and see if a unix executable file is generated in the current directory
  5. If yes, run the executable file using the following command, where 'cell_outline.csv' is a file with a single outline:
      ./loco-efa_example_GS cell_outline.csv
  6. In case of multiple files, store the outline files in a separate folder, for example Outlines, and loop over these files using the following line in the command line:
    for FILE in ./Outlines/*; do ./loco-efa_example_GS $FILE; done

Define which gastruloids are unique and plot their L2/L1 vs L3/L1 values (96h, 72h)¶

Place the numbers belonging to unique gastruloids in one sequence.¶

In [4]:
# Note: some of our images were overlapping each other, which caused some gastruloids to appear twice in our data set.
# We carefully checked all the extracted outlines and removed the doubles. Therefore, 500 - 413 = 87 extracted outlines were removed.

# 96 h gastruloids
a1 = np.arange(1,7)
a2 = np.arange(10,27)
a3 = np.arange(29,37)
a4 = np.array([38, 40])
a5 = np.arange(42,90)
a6 = np.arange(91,100)
a62 = np.array([102,104,106,107,109,112,118,121])
a7 = np.arange(123,157)
a8 = np.array([158,166])
a9 = np.arange(168,172)
a10 = np.array([173])
a11 = np.arange(175,183)
a12 = np.array([185,186,188])
a13 = np.arange(192,214)
a14 = np.arange(216,224)
a15 = np.array([225,226])
a16 = np.arange(228,243)
a17 = np.arange(244,252)
a18 = np.array([253,254])
a19 = np.arange(256,273)
a20 = np.array([274])
a21 = np.arange(276,282)
a22 = np.array([283])
a23 = np.arange(285,296)
a24 = np.array([297])
a25 = np.arange(299,304)
a26 = np.arange(306,311)
a27 = np.array([312,317,319,322,323])
a28 = np.arange(326,349)
a29 = np.arange(350,363)
a30 = np.arange(364,372)
a31 = np.arange(373,376)
a32 = np.arange(377,380)
a33 = np.arange(381,390)
a34 = np.arange(391,400)
a35 = np.arange(401,408)
a36 = np.arange(409,412)
a37 = np.arange(413,435)
a38 = np.arange(437,448)
a39 = np.arange(449,452)
a40 = np.arange(453,459)
a41 = np.arange(461,472)
a42 = np.array([473,474])
a43 = np.arange(476,485)
a44 = np.arange(486,490)
a45 = np.arange(492,500)

a_96h = np.concatenate((a1,a2,a3,a4,a5,a6,a62,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20,a21,a22,a23,a24,a25,a26,a27,a28,a29,a30,a31,a32,a33,a34,a35,a36,a37,a38,a39,a40,a41,a42,a43,a44,a45),dtype = np.int64)
print('Total number of 96 h gastruloids:', len(a_96h))

# 72 h gastruloids
a1 = np.arange(1,11)
a2 = np.arange(13,27)
a3 = np.array([30])
a4 = np.arange(32,89)
a5 = np.arange(91,101)
a6 = np.array([112])
a7 = np.array([114])
a8 = np.arange(116,127)
a9 = np.arange(130,134)
a10 = np.arange(136,140)
a11 = np.arange(146,147)
a12 = np.arange(149,155)
a13 = np.arange(159,171)

a_72h = np.concatenate((a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13))#,dtype = np.int64)
print('Total number of 72 h gastruloids:', len(a_72h))

# 10 uM ROCK inhibitor gastruloids
a1 = np.arange(1,19)
a2 = np.array([20])
a3 = np.arange(23,37)
a4 = np.array([38])
a5 = np.array([40])
a6 = np.arange(42,73)
a7 = np.arange(74,76)
a8 = np.arange(77,112)
a9 = np.array([113])
a10 = np.arange(115,118)
a11 = np.array([119])
a12 = np.array([121])
a13 = np.array([123])

a_10uM_RI = np.concatenate((a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13))#,dtype = np.int64)
print('ROCK inhibitor experiment - total number of treated gastruloids:', len(a_10uM_RI))

# No ROCK inhibitor gastruloids
a1 = np.arange(1,30)
a2 = np.arange(31,34)
a3 = np.arange(35,50)
a4 = np.arange(53,58)
a5 = np.arange(60,67)
a6 = np.arange(68,80)
a7 = np.arange(81,96)
a8 = np.arange(97,102)
a9 = np.arange(103,105)
a10 = np.arange(106,113)
a11 = np.array([115])
a12 = np.arange(117,122)

a_no_RI = np.concatenate((a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12))#,dtype = np.int64)
print('ROCK inhibitor experiment - total number of control gastruloids:', len(a_no_RI))
Total number of 96 h gastruloids: 413
Total number of 72 h gastruloids: 132
ROCK inhibitor experiment - total number of treated gastruloids: 110
ROCK inhibitor experiment - total number of control gastruloids: 106

Read in the txt file with L values for each gastruloid and divide L3 and L2 by L1 (normalization)¶

In [6]:
folder_72h = "/Users/esmeeadegeest/Documents/Outlines_72h_143unique/Outlines"
L2_div_L1_72h_array=[]
L3_div_L1_72h_array=[]
for i in range(np.size(a_72h)):
    i_string = str(a_72h[i])
    data_file = "Contour_"+i_string+"outline.txt"
    path_i = os.path.join(folder_72h,data_file)
    contents = pd.read_table(path_i)
    
    L_values_72h = pd.DataFrame(contents).to_numpy() # convert dataframe to array in numpy
    L2_div_L1_72h = np.divide(L_values_72h[1], L_values_72h[0]) # divide L2 by L1
    L3_div_L1_72h = np.divide(L_values_72h[2], L_values_72h[0]) # divide L3 by L1
    
    L2_div_L1_72h_array.append(L2_div_L1_72h) # store L2/L1 values
    L3_div_L1_72h_array.append(L3_div_L1_72h) # store L3/L1 values
In [7]:
folder_96h = "/Users/esmeeadegeest/Documents/L_values_no_doubles_413/"
L2_div_L1_96h_array=[]
L3_div_L1_96h_array=[]
for i in range(np.size(a_96h)):
    i_string = str(a_96h[i])
    data_file = "Data_"+i_string+".txt"
    path_i = os.path.join(folder_96h,data_file)
    contents = pd.read_table(path_i)
    
    L_values_96h = pd.DataFrame(contents).to_numpy() # convert dataframe to array in numpy
    L2_div_L1_96h = np.divide(L_values_96h[1], L_values_96h[0]) # divide L2 by L1
    L3_div_L1_96h = np.divide(L_values_96h[2], L_values_96h[0]) # divide L3 by L1
    
    L2_div_L1_96h_array.append(L2_div_L1_96h) # store L2/L1 values
    L3_div_L1_96h_array.append(L3_div_L1_96h) # store L3/L1 values
In [8]:
folder_No_RI = "/Users/esmeeadegeest/Documents/L-values_No-RI/"
L2_div_L1_No_RI_array=[]
L3_div_L1_No_RI_array=[]
for i in range(np.size(a_no_RI)):
    i_string = str(a_no_RI[i])
    data_file = "Contour_"+i_string+"outline.txt"
    path_i = os.path.join(folder_No_RI,data_file)
    contents = pd.read_table(path_i)
    
    L_values_No_RI = pd.DataFrame(contents).to_numpy() # convert dataframe to array in numpy
    L2_div_L1_No_RI = np.divide(L_values_No_RI[1], L_values_No_RI[0]) # divide L2 by L1
    L3_div_L1_No_RI = np.divide(L_values_No_RI[2], L_values_No_RI[0]) # divide L3 by L1
    
    L2_div_L1_No_RI_array.append(L2_div_L1_No_RI) # store L2/L1 values
    L3_div_L1_No_RI_array.append(L3_div_L1_No_RI) # store L3/L1 values
In [9]:
folder_10uM_RI = "/Users/esmeeadegeest/Documents/L-values_10uM-RI/"
L2_div_L1_10uM_RI_array=[]
L3_div_L1_10uM_RI_array=[]
for i in range(np.size(a_10uM_RI)):
    i_string = str(a_10uM_RI[i])
    data_file = "Contour_"+i_string+"outline.txt"
    path_i = os.path.join(folder_10uM_RI,data_file)
    contents = pd.read_table(path_i)
    
    L_values_10uM_RI = pd.DataFrame(contents).to_numpy() # convert dataframe to array in numpy
    L2_div_L1_10uM_RI = np.divide(L_values_10uM_RI[1], L_values_10uM_RI[0]) # divide L2 by L1
    L3_div_L1_10uM_RI = np.divide(L_values_10uM_RI[2], L_values_10uM_RI[0]) # divide L3 by L1
    
    L2_div_L1_10uM_RI_array.append(L2_div_L1_10uM_RI) # store L2/L1 values
    L3_div_L1_10uM_RI_array.append(L3_div_L1_10uM_RI) # store L3/L1 values

Plot the figures¶

In [10]:
#72 h and 96 h
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 
cbmap=['#377eb8', '#ff7f00', '#4daf4a',
                  '#f781bf', '#a65628', '#984ea3',
                  '#999999', '#e41a1c', '#dede00']

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

ax.scatter(L2_div_L1_72h_array,L3_div_L1_72h_array, s=1, marker ="o", label='72h gastruloids', c=cbmap[5])
ax.scatter(L2_div_L1_96h_array,L3_div_L1_96h_array, s=1, marker ="o", label='96h gastruloids', c=cbmap[1])

fig.set_size_inches(10/2.54, 8/2.54)
plt.xlabel('$L_2/L_1$')
plt.ylabel('$L_3/L_1$')
#plt.title('all <-> all')
handles, labels = ax.get_legend_handles_labels()
lgd = ax.legend(handles, labels, loc='upper left', bbox_to_anchor=(1,1))
fig.set_size_inches(10/2.54, 8/2.54)

# save figure
#fig.savefig('Figures_paper/Figure_1_LOCO-EFA_72h_96h_purple_orange.svg', bbox_extra_artists=([lgd]), bbox_inches='tight')
In [11]:
# ROCK inhibitor experiment
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 
cbmap=['#377eb8', '#ff7f00', '#4daf4a',
                  '#f781bf', '#a65628', '#984ea3',
                  '#999999', '#e41a1c', '#dede00']

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

ax.scatter(L2_div_L1_No_RI_array,L3_div_L1_No_RI_array, s=1, marker ="o", label='Control', c=cbmap[1])
ax.scatter(L2_div_L1_10uM_RI_array,L3_div_L1_10uM_RI_array, s=1, marker ="o", label=r'10 $\mu$M Y-27632', c=cbmap[5])

fig.set_size_inches(10/2.54, 8/2.54)
plt.xlabel('$L_2/L_1$')
plt.ylabel('$L_3/L_1$')
plt.xlim(-0.02, 0.7)
plt.ylim(-0.015, 0.305)
#plt.title('all <-> all')
handles, labels = ax.get_legend_handles_labels()
lgd = ax.legend(handles, labels, loc='upper left', bbox_to_anchor=(1,1))
fig.set_size_inches(10/2.54, 8/2.54)

# save figure
#fig.savefig('Documents/Figure_X_LOCO-EFA_RI_purple_orange.svg', bbox_extra_artists=([lgd]), bbox_inches='tight')

In vitro time lapse analysis¶

0a. Create binary image¶

Our time lapse data consists of 31 image stacks, corresponding to 31 time points. Here, we used maximum projections in Z to analyze the 3D data in 2D. To determine where cells are with respect to the long and short axis of the gastruloid, we fit an ellipse to the 2D gastruloid shape. To obtain the 2D gastruloid shape, we first pre-process the images in FIJI to enhance the contrast. Then we use several morphological operations to create a binary image that shows the shape of the gastruloid. The binary image is a rough estimate of the gastruloid's shape, since not all cells within the gastruloid were mCherry+. From the binary images, an outline is created, which is used to fit an ellipse to.

In [3]:
# Pre-processing in FIJI:
# run("RGB Color");
# run("Enhance Contrast", "saturated=0.35");
# should be stored as PNG file.)

%matplotlib inline
img = []
img_GB = []
img_GB_thresh = []
img_binary = []

plt.rcParams['figure.figsize'] = [30, 20]
kernel = np.ones((5, 5), np.uint8) # Taking a matrix of size 5 as the kernel for the morphological operations

# Perform morphological operations on the maximum z projection for each time point using a specified threshold value and store the images in img_binary
for i in range(1,32):
    img_read = cv2.imread('Mosaic_TL_mCherry/Max_projections_Z/MAX_T'+str(i)+'_TL-P1-10min-33ts_G2_561nm-1700mA_160ms_6p7_4um_40zimage1_crop.png',0)
    img.append(img_read)
    
    if i == 1:
        threshold_value = 63 # We used 40 for all frames, except for frame of t = 1 (63) and t = 26 (45).
    elif i == 26:
        threshold_value = 45
    else:
        threshold_value = 40
        
    # Apply Gaussian blur
    img_read_GB = cv2.GaussianBlur(img_read, (51,51), cv2.BORDER_DEFAULT) # kernel size should be odd
    img_GB.append(img_read_GB)
    # Apply threshold value >> check in the next cell for each time frame whether the threshold value is correct
    ret,img_GB_read_thresh1 = cv2.threshold(img_read_GB,threshold_value,255,cv2.THRESH_BINARY)
    img_GB_thresh.append(img_GB_read_thresh1)
    # Apply dilation multiple times to fill the holes. 
    img_dilation = cv2.dilate(img_GB_read_thresh1, kernel, iterations=40)
    # Apply subsequently erosion to retrieve the initial size
    img_erosion = cv2.erode(img_dilation, kernel, iterations=40)
    img_binary.append(img_erosion)
In [4]:
%matplotlib inline
plt.rcParams['figure.figsize'] = [15, 15]
fig, axs = plt.subplots(2, 2)

# Show for one time point the original image, the image after Gaussian blur, 
# the image after applying a threshold, and the image after dilation and erosion.
n = 31 # time point
axs[0,0].imshow(img[n-1], cmap='gray')
axs[0,1].imshow(img_GB[n-1], cmap='gray')
axs[1,0].imshow(img_GB_thresh[n-1], cmap='gray')
axs[1,1].imshow(img_binary[n-1], cmap='gray')

plt.show()
In [ ]:
# Save binary images, note down the threshold value for each image
# We used threshold value 40 for all frames, except for the frame of t = 1 (63) and t = 26 (45).

# t = 1, threshold = 63
cv2.imwrite('Mosaic_TL_mCherry/Binary_images/t1-63-binary.png', img_binary[0])

# t = 2-25, threshold = 40
for i in range(2,26):
    cv2.imwrite('Mosaic_TL_mCherry/Binary_images/t'+str(i)+'-40-binary.png', img_binary[i-1])
    
# t = 26, threshold = 45
cv2.imwrite('Mosaic_TL_mCherry/Binary_images/t26-45-binary.png', img_binary[25])

# t = 27-31, threshold = 40
for i in range(27,32):
    cv2.imwrite('Mosaic_TL_mCherry/Binary_images/t'+str(i)+'-40-binary.png', img_binary[i-1])
In [5]:
# Area calculation and comparison t1 and t31
image_loc_t1 = 'Mosaic_TL_mCherry/Binary_images/t1-63-binary.png'
image_loc_t31 = 'Mosaic_TL_mCherry/Binary_images/t31-40-binary.png'

mask_t1 = cv2.imread(image_loc_t1)
mask_t31 = cv2.imread(image_loc_t31)

imgray_t1 = cv2.cvtColor(mask_t1,cv2.COLOR_BGR2GRAY)
imgray_t31 = cv2.cvtColor(mask_t31,cv2.COLOR_BGR2GRAY)

area_t1 = np.sum(imgray_t1)/255
area_t31 = np.sum(imgray_t31)/255
print('Ratio area t31/t1:', area_t31/area_t1)

plt.rcParams['figure.figsize'] = [10, 10]
fig, axs = plt.subplots(1, 2)
axs[0].imshow(mask_t1, cmap='gray')
axs[1].imshow(mask_t31, cmap='gray')
Ratio area t31/t1: 1.3836912100209595
Out[5]:
<matplotlib.image.AxesImage at 0x7fb20addab80>

0b. Extract outlines from binary images¶

In [ ]:
# We used the following code to extract outlines from the binary images.
def extract_outline_iv():

    for j in range(1,32):
        ind = list() # index
        cont_all = list() # contours
        
        if j == 1:
            image_loc = 'Mosaic_TL_mCherry/Binary_images/t'+str(j)+'-63-binary.png'
            file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-63-binary_outline.txt'
        elif j == 26:
            image_loc = 'Mosaic_TL_mCherry/Binary_images/t'+str(j)+'-45-binary.png'
            file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-45-binary_outline.txt'
        else:
            image_loc = 'Mosaic_TL_mCherry/Binary_images/t'+str(j)+'-40-binary.png'
            file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-40-binary_outline.txt'
            
        file1 = open(file_loc,"a")
        mask = cv2.imread(image_loc)

        if mask is not None:

            imgray = cv2.cvtColor(mask,cv2.COLOR_BGR2GRAY)
            ret,thresh = cv2.threshold(imgray,200,255,0)
            contours, im2 = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_NONE)         
            maxsizeloc = 0
            maxsize = 0
           
            for i in range (len(contours)):
                if contours[i].shape[0] > maxsize:
                    maxsizeloc = i
                    maxsize = contours[i].shape[0]
            cont = np.zeros((contours[maxsizeloc].shape[0], 2))
            cont[:,0] = contours[maxsizeloc][:,0,0]
            cont[:,1] = contours[maxsizeloc][:,0,1]

            plt.plot(cont[:,0], cont[:,1])
            cont_all.append(cont)
            
            # Write outlines in txt file
            for i in range (maxsize):
                file1.write(str(cont_all[0][i][0]) + ',' + str(cont_all[0][i][1]) + '\n')
            file1.close()
          
    return(cont_all)

extract_outline_iv()

1. Time lapse: load outlines, fit ellipse and store ellipse parameters¶

In [6]:
# Ellipse fit parameters
xc_iv = np.empty(31) # center coordinate x of ellipse
yc_iv = np.empty(31) # center coordinate y of ellipse
a_iv = np.empty(31) # length short arm of ellipse
b_iv = np.empty(31) # length long arm of ellipse
theta_iv= np.empty(31) # rotation angle of ellipse

# Figure
%matplotlib inline
fig, axs = plt.subplots(2, 16,figsize=(100,30))#sharex=True, sharey=True)

# Load outline coordinates for each time step and fit ellipse, store ellipse parameters
for j in range (1,32):
    # txt file name in which outline is stored 
    if j == 1:
        outline_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-63-binary_outline.txt'
        file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-63_xycenter-theta-axes.txt' # assign a file location and open it to later on write the ellipse parameters in the file
    elif j == 26:
        outline_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-45-binary_outline.txt'
        file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-45_xycenter-theta-axes.txt'
    else:
        outline_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-40-binary_outline.txt'
        file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-40_xycenter-theta-axes.txt'
   
    file1 = open(file_loc,"a")
    
    # open outline file and store coordinates in array
    with open(outline_loc) as file_name:
        points = np.loadtxt(file_name, delimiter=",")
    a_points = np.array(points)
    x = a_points[:, 0]
    y = a_points[:, 1]
    
    # fit ellipse to array and get parameters
    ell = EllipseModel()
    ell.estimate(a_points)
    xc_iv[j-1], yc_iv[j-1], a_iv[j-1], b_iv[j-1], theta_iv[j-1] = ell.params
    
    # Write parameters to txt file
    #file1.write(str(xc_iv[j-1]) + ',' + str(yc_iv[j-1]) + ',' + str(theta_iv[j-1]) + ',' + str(a_iv[j-1]) + ',' + str(b_iv[j-1]) + '\n')
    
    # Plot ellipse fits
    if j < 16: # top row of images
        axs[0,j-1].scatter(x, y) # outline
        axs[0,j-1].scatter(xc_iv[j-1], yc_iv[j-1], color='red', s=100) # center of ellipse
        axs[0,j-1].set_xlim(0,1600)
        axs[0,j-1].set_ylim(0,3000)
        axs[0,j-1].set_xticks([])
        axs[0,j-1].set_yticks([])
        
        # ellipse
        ell_patch = Ellipse((xc_iv[j-1], yc_iv[j-1]), 2*a_iv[j-1], 2*b_iv[j-1], theta_iv[j-1]*180/np.pi, edgecolor='red', facecolor='none')
        # ellipse with orientation angle 0
        ell_patch2 = Ellipse((xc_iv[j-1], yc_iv[j-1]), 2*a_iv[j-1], 2*b_iv[j-1], 0, edgecolor='black', facecolor='none')

        axs[0,j-1].add_patch(ell_patch)
        axs[0,j-1].add_patch(ell_patch2)
        
    else: # bottom row of images
        axs[1,j-16].scatter(x, y) # outline
        axs[1,j-16].scatter(xc_iv[j-1], yc_iv[j-1], color='red', s=100) # center of ellipse
        axs[1,j-16].set_xlim(0,1600)
        axs[1,j-16].set_ylim(0,3000)
        axs[1,j-16].set_xticks([])
        axs[1,j-16].set_yticks([])
        
        # ellipse
        ell_patch = Ellipse((xc_iv[j-1], yc_iv[j-1]), 2*a_iv[j-1], 2*b_iv[j-1], theta_iv[j-1]*180/np.pi, edgecolor='red', facecolor='none')
        # ellipse with orientation angle 0
        ell_patch2 = Ellipse((xc_iv[j-1], yc_iv[j-1]), 2*a_iv[j-1], 2*b_iv[j-1], 0, edgecolor='black', facecolor='none')

        axs[1,j-16].add_patch(ell_patch)
        axs[1,j-16].add_patch(ell_patch2)
plt.show()

# From the images we can tell that for some shapes the short arm length and long arm length are
# stored in the wrong parameters: b and a, respectively, while they should be stored in a and b, respectively.
In [7]:
# Here we print the short arm length (a) and long arm length (b) and see that they are not always stored correctly.
print('Time point 1 (correct): a = ',a_iv[0],'b =', b_iv[0])
print('Time point 4 (incorrect): a = ',a_iv[3],'b = ', b_iv[3])
Time point 1 (correct): a =  588.7241640888858 b = 879.3560145828596
Time point 4 (incorrect): a =  921.3089122428711 b =  591.6155040634976

2. Cell tracking: convert ROIs to arrays and store in list¶

In [8]:
# We traced cells manually using the point tool in ImageJ. The trajectories were stored in an ROI set.
# Import the roi_set of each cell
N_cells = 14
cell_list = [None] * N_cells
roi_set = [None] * N_cells

for i in range(N_cells):
    path_file = 'Mosaic_TL_mCherry/Cells_by_number/RoiSet_Cell_' + str(i+1) + '.zip'
    roi_set[i] = read_roi_zip(path_file)

We also had dividing cells and therefore tracked the daughter cells after division. One cell was tracked starting from the last time point to the first time point for convenience. The following cells have different scripts to convert their ROIs to arrays:

cell 3 = dividing cell, positions of daughter cells alternate
cell 4 = dividing cell, position of first daughter cell first
cell 5 = dividing cell, position of first daughter cell first
cell 7 = dividing cell, position of first daughter cell first
cell 8 = order in set is from last to first time point
cell 9 = dividing cell, position of first daughter cell first

For each cell, we store the trajectories in a list. Depending on whether a cell divides, we store the trajectory of the mother cell in path_x, path_y and the trajectories of the daughter cells in path_x_d1, path_y_d1 and path_x_d2, path_y_d2.
cell_list = path_x, path_y, division_status, time_step_division, path_x_d1, path_y_d1, path_x_d2, path_y_d2

In [9]:
# Dividing cell script for cell 3 (alternating position)
n = 3
roi_set_sel = roi_set[n-1]
time_step_division = 17
t = 1 # start time

division_status = 'Dividing cell'
path_x = np.arange(0) # mother cell x
path_y = np.arange(0) # mother cell y
path_x_d1 = np.arange(0) # daughter cell 1 x
path_x_d2 = np.arange(0) # daughter cell 2 x
path_y_d1 = np.arange(0) # daughter cell 1 y
path_y_d2 = np.arange(0) # daughter cell 2 xy
roi_set_length = len(roi_set_sel)
d1 = bool(True) # daughter cell 1 is set to true

for key, value in roi_set_sel.items():
    path_roi = roi_set_sel[key]
    
    if t < time_step_division:
        path_x = np.append(path_x, path_roi['x'][0])
        path_y = np.append(path_y, path_roi['y'][0])
        t = t + 1
    elif t >= time_step_division:
        if d1 == True:
            path_x_d1 = np.append(path_x_d1, path_roi['x'][0])
            path_y_d1 = np.append(path_y_d1, path_roi['y'][0])
            d1 = bool(False) # set daughter cell 1 to false > store next position in array of 2nd daughter cell
            t = t + 1
        elif d1 == False:
            path_x_d2 = np.append(path_x_d2, path_roi['x'][0])
            path_y_d2 = np.append(path_y_d2, path_roi['y'][0])
            d1 = bool(True) # set daughter cell 1 to true > store next position in array of 1st daughter cell
            t = t + 1
            
# store the trajectories in a list   
cell_list[n-1] = [path_x, path_y, division_status, time_step_division, path_x_d1[:(len(path_x_d1))], path_y_d1[:(len(path_y_d1))], path_x_d2[:(len(path_x_d2))], path_y_d2[:(len(path_y_d2))]]
In [10]:
# Dividing cell script for cell 5 (weird order: first alternating positions, then paths subsequently stored in ROI set)
n = 5
roi_set_sel = roi_set[n-1]
time_step_division_par = 9
time_step_division_sep = 16
t_sep_start = time_step_division_par + 2*(time_step_division_sep - time_step_division_par)
t_sep_d2_start = t_sep_start + (31 - time_step_division_sep + 1)

roi_set_length = len(roi_set_sel)
t = 1
division_status = 'Dividing cell'
path_x = np.arange(0)
path_y = np.arange(0)
path_x_d1 = np.arange(0)
path_x_d2 = np.arange(0)
path_y_d1 = np.arange(0)
path_y_d2 = np.arange(0)
d1 = bool(True)

for key, value in roi_set_sel.items():
    path_roi = roi_set_sel[key]
    
    if t < time_step_division_par: # to 8 (incl)
        path_x = np.append(path_x, path_roi['x'][0])
        path_y = np.append(path_y, path_roi['y'][0])
        t = t + 1

    elif t >= time_step_division_par and t < t_sep_start: # 9 to 15 (incl)
        if d1 == True:
            path_x_d1 = np.append(path_x_d1, path_roi['x'][0])
            path_y_d1 = np.append(path_y_d1, path_roi['y'][0])
            d1 = bool(False)
            t = t + 1
        elif d1 == False:
            path_x_d2 = np.append(path_x_d2, path_roi['x'][0])
            path_y_d2 = np.append(path_y_d2, path_roi['y'][0])
            d1 = bool(True)
            t = t + 1
    elif t >= t_sep_start and t < t_sep_d2_start:
        path_x_d1 = np.append(path_x_d1, path_roi['x'][0])
        path_y_d1 = np.append(path_y_d1, path_roi['y'][0])
        t = t + 1
    elif t >= t_sep_d2_start:
        path_x_d2 = np.append(path_x_d2, path_roi['x'][0])
        path_y_d2 = np.append(path_y_d2, path_roi['y'][0])
        t = t + 1

cell_list[n-1] = [path_x, path_y, division_status, time_step_division_par, path_x_d1[:(len(path_x_d1))], path_y_d1[:(len(path_y_d1))], path_x_d2[:(len(path_x_d2))], path_y_d2[:(len(path_y_d2))]]
In [11]:
# Dividing cell script for cell 7
n = 7
time_step_division = 14
roi_set_sel = roi_set[n-1]
roi_set_length = len(roi_set_sel)

t = 1
division_status = 'Dividing cell'
path_x = np.arange(0)
path_y = np.arange(0)
path_x_d1 = np.arange(0)
path_x_d2 = np.arange(0)
path_y_d1 = np.arange(0)
path_y_d2 = np.arange(0)

for key, value in roi_set_sel.items():
    path_roi = roi_set_sel[key]
    if t < time_step_division:
        path_x = np.append(path_x, path_roi['x'][0])
        path_y = np.append(path_y, path_roi['y'][0])
        t = t + 1
    elif t >= time_step_division and t <= 31:
        path_x_d1 = np.append(path_x_d1, path_roi['x'][0])
        path_y_d1 = np.append(path_y_d1, path_roi['y'][0])
        t = t + 1
    elif t > 31:
        path_x_d2 = np.append(path_x_d2, path_roi['x'][0])
        path_y_d2 = np.append(path_y_d2, path_roi['y'][0])
        t = t + 1

cell_list[n-1] = [path_x, path_y, division_status, time_step_division, path_x_d1[:(len(path_x_d1))], path_y_d1[:(len(path_y_d1))], path_x_d2[:(len(path_x_d2))], path_y_d2[:(len(path_y_d2))]]
In [12]:
# Dividing cell script for cell 9
n = 9
time_step_division = 12
roi_set_sel = roi_set[n-1]

roi_set_length = len(roi_set_sel)
t = 1
division_status = 'Dividing cell'
path_x = np.arange(0)
path_y = np.arange(0)
path_x_d1 = np.arange(0)
path_x_d2 = np.arange(0)
path_y_d1 = np.arange(0)
path_y_d2 = np.arange(0)

for key, value in roi_set_sel.items():
    path_roi = roi_set_sel[key]
    
    if t < time_step_division:
        path_x = np.append(path_x, path_roi['x'][0])
        path_y = np.append(path_y, path_roi['y'][0])
        t = t + 1

    elif t >= time_step_division and t <= 31:
        path_x_d1 = np.append(path_x_d1, path_roi['x'][0])
        path_y_d1 = np.append(path_y_d1, path_roi['y'][0])
        t = t + 1

    elif t > 31:
        path_x_d2 = np.append(path_x_d2, path_roi['x'][0])
        path_y_d2 = np.append(path_y_d2, path_roi['y'][0])
        t = t + 1

cell_list[n-1] = [path_x, path_y, division_status, time_step_division, path_x_d1[:(len(path_x_d1))], path_y_d1[:(len(path_y_d1))], path_x_d2[:(len(path_x_d2))], path_y_d2[:(len(path_y_d2))]]
In [13]:
# Non-dividing cells: 1, 2, 4, 6, 10, 11, 12, 13, 14
# Fill in here:
n = [1, 2, 4, 6, 10, 11, 12, 13, 14]

division_status = 'Non-dividing cell'

for i in n:
    path_x = np.arange(0)
    path_y = np.arange(0)
    roi_set_sel = roi_set[i-1]
    
    for key, value in roi_set_sel.items():
        path_roi = roi_set_sel[key]
        path_x = np.append(path_x,path_roi['x'][0])
        path_y = np.append(path_y,path_roi['y'][0])
        
    cell_list[i-1] = [path_x[:(len(path_x))], path_y[:(len(path_y))], division_status]
In [14]:
# Non-dividing cell - cell 8 (reverse order)
# Fill in here:
n = 8
roi_set_sel = roi_set[n-1]

division_status = 'Non-dividing cell'
path_x = np.arange(0)
path_y = np.arange(0)

for key, value in roi_set_sel.items():
    path_roi = roi_set_sel[key]
    path_x = np.append(path_x,path_roi['x'][0])
    path_y = np.append(path_y,path_roi['y'][0])

path_x_rev = path_x[::-1]
path_y_rev = path_y[::-1]

cell_list[n-1] = [path_x_rev[:(len(path_x_rev))], path_y_rev[:(len(path_y_rev))], division_status]

3. Analysis: rotate the positions of the cells to calculate their position w.r.t. the long axis¶

In [15]:
n_cells = 14 
divcells = 4 # 4 dividing cells
#cells_sel = [1, 2, 4, 6, 8, 10, 11, 12, 13, 14] #non-dividing cells only
cells_sel = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14] #all cells
div_cells_no = [3, 5, 7, 9]

t_steps = 31
go_twice = False # is true for dividing cells
count = 0
theta_corr_iv = [None] * t_steps

rotated_cell_centers_iv = np.empty((n_cells+divcells, t_steps, 4)) # we store: x, y, V_x, V_y

for i in range(n_cells):

    sel_cell = cells_sel[i]-1
        
    if cell_list[sel_cell][2] == 'Dividing cell':
       
        go_twice = True
    
        x_d1 = []
        y_d1 = []
        x_d2 = []
        y_d2 = []

        x_d1 = np.append([cell_list[sel_cell][0]],[cell_list[sel_cell][4]]) # paste trajectory (x coordinate) daughter cell 1 after trajectory mother cell
        y_d1 = np.append([cell_list[sel_cell][1]],[cell_list[sel_cell][5]]) # paste trajectory (y coordinate) daughter cell 1 after trajectory mother cell
    
        x_d2 = np.append([cell_list[sel_cell][0]],[cell_list[sel_cell][6]]) # paste trajectory (x coordinate) daughter cell 2 to trajectory mother cell
        y_d2 = np.append([cell_list[sel_cell][1]],[cell_list[sel_cell][7]]) # paste trajectory (y coordinate) daughter cell 2 to trajectory mother cell
    
    else:
        go_twice = False
    
    if go_twice == False: # non-dividing cell

        for j in range(t_steps):

            if theta_iv[j] > 1.5708: # if gastruloid (ellipse fit) orientation angle is larger than pi/2, correct the angle by subtracting pi/2
                theta_corr_iv[j] = theta_iv[j] - 1.5708
            else:
                theta_corr_iv[j] = theta_iv[j]
            
            # 1. Calculate the difference in x and y position of a cell w.r.t. the center of the ellipse for each time step
            dx = cell_list[sel_cell][0][j]-xc_iv[j]
            dy = cell_list[sel_cell][1][j]-yc_iv[j]
            
            # 2. Calculate the difference in x and y position of each cell w.r.t. to the new x and y axes, the axes of the ellipse, where the long axis is 
            # the new y axis and the short axis is the new x axis. We call them V_x and V_y.
            V_x = np.cos(theta_corr_iv[j]) * dx + np.sin(theta_corr_iv[j]) * dy
            V_y = np.cos(theta_corr_iv[j]) * dy - np.sin(theta_corr_iv[j]) * dx
            
            # 3. Store V_x and V_y and store the new x and y coordinates: V_x + x(center of ellipse) and V_y + y(center of ellipse)
            rotated_cell_centers_iv[i+count,j,0] = V_x + xc_iv[j]
            rotated_cell_centers_iv[i+count,j,2] = V_x
            
            rotated_cell_centers_iv[i+count,j,1] = V_y + yc_iv[j]
            rotated_cell_centers_iv[i+count,j,3] = V_y
            
                
    elif go_twice == True: # dividing cell
        
        for j in range(t_steps):
            
            if theta_iv[j] > 1.5708: # if gastruloid (ellipse fit) orientation angle is larger than pi/2, correct the angle by subtracting pi/2
                theta_corr_iv[j] = theta_iv[j] - 1.5708
            else:
                theta_corr_iv[j] = theta_iv[j]
              
            # 1. Calculate the difference in x and y position of a cell w.r.t. the center of the ellipse for each time step
            dx = x_d1[j]-xc_iv[j]
            dy = y_d1[j]-yc_iv[j]
            
            # 2. Calculate the difference in x and y position of each cell w.r.t. to the new x and y axes, the axes of the ellipse, where the long axis is 
            # the new y axis and the short axis is the new x axis. We call them V_x and V_y.
            V_x = np.cos(theta_corr_iv[j]) * dx + np.sin(theta_corr_iv[j]) * dy
            V_y = np.cos(theta_corr_iv[j]) * dy - np.sin(theta_corr_iv[j]) * dx
            
            # 3. Store V_x and V_y and store the new x and y coordinates: V_x + x(center of ellipse) and V_y + y(center of ellipse)
            rotated_cell_centers_iv[i+count,j,0] = V_x + xc_iv[j]
            rotated_cell_centers_iv[i+count,j,2] = V_x
            
            rotated_cell_centers_iv[i+count,j,1] = V_y + yc_iv[j]
            rotated_cell_centers_iv[i+count,j,3] = V_y
        
        count = count + 1 # for 2nd daughter cell
        
        for j in range(t_steps):
            
            if theta_iv[j] > 1.5708: # if gastruloid (ellipse fit) orientation angle is larger than pi/2, correct the angle by subtracting pi/2
                theta_corr_iv[j] = theta_iv[j] - 1.5708
            else:
                theta_corr_iv[j] = theta_iv[j]
            
            # 1. Calculate the difference in x and y position of a cell w.r.t. the center of the ellipse for each time step
            dx = x_d2[j]-xc_iv[j]
            dy = y_d2[j]-yc_iv[j]
            
            # 2. Calculate the difference in x and y position of each cell w.r.t. to the new x and y axes, the axes of the ellipse, where the long axis is 
            # the new y axis and the short axis is the new x axis. We call them V_x and V_y.
            V_x = np.cos(theta_corr_iv[j]) * dx + np.sin(theta_corr_iv[j]) * dy
            V_y = np.cos(theta_corr_iv[j]) * dy - np.sin(theta_corr_iv[j]) * dx
            
            # 3. Store V_x and V_y and store the new x and y coordinates: V_x + x(center of ellipse) and V_y + y(center of ellipse)
            rotated_cell_centers_iv[i+count,j,0] = V_x + xc_iv[j]
            rotated_cell_centers_iv[i+count,j,2] = V_x
            
            rotated_cell_centers_iv[i+count,j,1] = V_y + yc_iv[j]
            rotated_cell_centers_iv[i+count,j,3] = V_y    

4. Plot: starting points of 14 cells on t = 1 image¶

In [16]:
# T = 1, 14 cells starting points

%matplotlib inline
plt.rcParams['figure.figsize'] = [30, 20]
imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t1-200epochsRGB_SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])

# Colors
cmap = [#'#000000', #black
        '#e69d00', #orange
        '#f0e442', #yellow
        '#009e73', #bluish green
        '#56b4e9', #sky blue
        '#cc79a7', #reddish purple
        '#d55e00', #vermillion (reddish)
        '#0072b2', #blue
        #'#ffffff', #white
        '#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$']
n = 14

for i in range(n):
    if i >= 9:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
       
    i_sel = i
    cell_marker_i = cell_marker[i_sel]
    if i < 8:
        color_i = cmap[i]
    if i >= 8 and i < 16:
        color_i = cmap[i-8] 

    # Plot starting points 
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker='o',
     linewidth=2, markersize=23) # start point black edge
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = color_i, marker='o',
     linewidth=2, markersize=21) # start point colored
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'w', marker='o',
     linewidth=2, markersize=16) # start point white
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker=cell_marker_i,
     linewidth=2, markersize=10*multiplying_factor) # start point number
          
#get current axes
ax = plt.gca()

#hide x-axis
ax.get_xaxis().set_visible(False)

#hide y-axis 
ax.get_yaxis().set_visible(False)

# Save figure
#plt.savefig("Mosaic_TL_mCherry/Figure_in_vitro_tracking_mosaic_icons_t1_more_white.png", dpi=600, bbox_inches='tight',pad_inches = 0)

# Plot figure
fig = plt.show()

5a. Plot: trajectories of 14 cells on t = 31 image¶

In [17]:
# T = 31, 14 cells trajectories
%matplotlib inline
plt.rcParams['figure.figsize'] = [30, 20]

imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t31-200epochsRGB-SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])

# Colors
cmap = [#'#000000', #black
        '#e69d00', #orange
        '#f0e442', #yellow
        '#009e73', #bluish green
        '#56b4e9', #sky blue
        '#cc79a7', #reddish purple
        '#d55e00', #vermillion (reddish)
        '#0072b2', #blue
        #'#ffffff', #white
        '#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$']
n = 14

for i in range(n):
    if i >= 9:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
       
    i_sel = i
    cell_marker_i = cell_marker[i_sel]
    if i < 8:
        color_i = cmap[i]
    if i >= 8 and i < 16:
        color_i = cmap[i-8]
    
    # dividing cells
    # division points = cell 3: 17, cell 5: 16, cell 7: 14, cell 9: 12
    if cell_list[i_sel][2] == 'Dividing cell':
        cell_division_point = cell_list[i_sel][3]
        
        plt.plot(cell_list[i_sel][4], cell_list[i_sel][5], color = 'k', #linestyle = 'dotted',
         linewidth = 3, markersize = 3) # black trajectory "daughter cell 1"
        plt.plot(cell_list[i_sel][4], cell_list[i_sel][5], color = color_i,
          linewidth=2.5, markersize=2.5) # colored trajectory "daughter cell 1"
        
        plt.plot(cell_list[i_sel][6], cell_list[i_sel][7], color = 'k', #linestyle = 'dotted',
         linewidth = 3, markersize = 3) # black trajectory "daughter cell 2"
        plt.plot(cell_list[i_sel][6], cell_list[i_sel][7], color = color_i,
          linewidth=2.5, markersize=2.5) # colored trajectory "daughter cell 2"
    
    # Trajectory from start to end (mother cell or non-dividing cell)
    plt.plot(cell_list[i_sel][0], cell_list[i_sel][1], color = 'k', #linestyle = 'dotted',
     linewidth = 3, markersize = 3) # line connecting start and end point black 
    plt.plot(cell_list[i_sel][0], cell_list[i_sel][1], color = color_i,
     linewidth=2.5, markersize=2.5) # line connecting start and end point color

    # start points
    if cell_list[i_sel][2] == 'Dividing cell':
        
        # x and y coordinates for connecting line mother cell and daughter cell
        mother_x = cell_list[i_sel][0][-1]
        mother_y = cell_list[i_sel][1][-1]
        daughter_1_x = cell_list[i_sel][4][0]
        daughter_1_y = cell_list[i_sel][5][0]
        daughter_2_x = cell_list[i_sel][6][0]
        daughter_2_y = cell_list[i_sel][7][0]
        x_points_d1 = [mother_x, daughter_1_x]
        x_points_d2 = [mother_x, daughter_2_x]
        y_points_d1 = [mother_y, daughter_1_y]
        y_points_d2 = [mother_y, daughter_2_y]
        
        # plot connecting line for daughter cell 1
        plt.plot(x_points_d1, y_points_d1, color = 'k',# linestyle = 'dotted',
         linewidth = 3, markersize = 3) # line connecting start and end point black 
        plt.plot(x_points_d1, y_points_d1, color = color_i,
         linewidth=2.5, markersize=2.5) # line connecting start and end point color
        
        # plot connecting line for daughter cell 2
        plt.plot(x_points_d2, y_points_d2, color = 'k', #linestyle = 'dotted',
         linewidth = 3, markersize = 3) # line connecting start and end point black 
        plt.plot(x_points_d2, y_points_d2, color = color_i,
         linewidth=2.5, markersize=2.5) # line connecting start and end point color
        
        # start point daughter cell 1
        plt.plot(cell_list[i_sel][4][0],cell_list[i_sel][5][0], color = 'k',
          marker='o', linewidth=20, markersize=23) # end point black edge
        plt.plot(cell_list[i_sel][4][0],cell_list[i_sel][5][0], color = color_i,
          marker='o', linewidth=20, markersize=21) # end point colored
        plt.plot(cell_list[i_sel][4][0],cell_list[i_sel][5][0], color = 'k',
          marker='o', linewidth=20, markersize=16) #14 end point black
        plt.plot(cell_list[i_sel][4][0],cell_list[i_sel][5][0], color = 'w', marker=cell_marker_i,
         linewidth=10, markersize=10*multiplying_factor) # end point number in white
        
        # start point daughter cell 2
        plt.plot(cell_list[i_sel][6][0],cell_list[i_sel][7][0], color = 'k',
          marker='o', linewidth=20, markersize=23) # end point black edge
        plt.plot(cell_list[i_sel][6][0],cell_list[i_sel][7][0], color = color_i,
          marker='o', linewidth=20, markersize=21) # end point colored
        plt.plot(cell_list[i_sel][6][0],cell_list[i_sel][7][0], color = 'k',
          marker='o', linewidth=20, markersize=16) # end point black
        plt.plot(cell_list[i_sel][6][0],cell_list[i_sel][7][0], color = 'w', marker=cell_marker_i,
         linewidth=20, markersize=10*multiplying_factor) # end point number
    
    # start points non-dividing cells and mother cells
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker='o',
     linewidth=2, markersize=23) # start point black edge
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = color_i, marker='o',
     linewidth=2, markersize=21) # start point colored
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'w', marker='o',
     linewidth=2, markersize=16) # start point white
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker=cell_marker_i,
     linewidth=2, markersize=10*multiplying_factor) # start point number
    
    # end points
    end_points = 1
    if end_points == 1:
        if cell_list[i_sel][2] == 'Dividing cell':
            # daughter cell 1
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) # end point white
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=8*multiplying_factor) # end point number
        
            # daughter cell 2
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) #13 end point white
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=10*multiplying_factor) # end point number
        
        else: # no dividing cell
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) # end point white
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=10*multiplying_factor) # end point number

#get current axes
ax = plt.gca()

#hide x-axis
ax.get_xaxis().set_visible(False)

#hide y-axis 
ax.get_yaxis().set_visible(False)

# Save figure
#plt.savefig("Mosaic_TL_mCherry/Figure_in_vitro_tracking_mosaic_icons-t31_continuous-line_black.png", dpi=600, bbox_inches='tight',pad_inches = 0)

# Plot figure
fig = plt.show()

5b. Plot: correct trajectory plot for drift of gastruloid only¶

In [18]:
# CORRECTION FOR DRIFT only
# For T = 31, plot 14 cell trajectories that are corrected for the 'drift' of the center's gastruloid (center of ellipse) over time
# With drift we mean: how the center of the gastruloid - estimated by the center of the fitted ellipse - changes in position over time w.r.t. to the center of the gastruloid at the final time step.
%matplotlib inline
plt.rcParams['figure.figsize'] = [30, 20]

imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t31-200epochsRGB-SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])

cmap = [#'#000000', #black
        '#e69d00', #orange
        '#f0e442', #yellow
        '#009e73', #bluish green
        '#56b4e9', #sky blue
        '#cc79a7', #reddish purple
        '#d55e00', #vermillion (reddish)
        '#0072b2', #blue
        #'#ffffff', #white
        '#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$']
n = 14

# Calculate drift
x_corr_drift_iv = [None] * t_steps
y_corr_drift_iv = [None] * t_steps
for j in range(0,t_steps):
    x_corr_drift_iv[j] = xc_iv[-1]-xc_iv[j]
    y_corr_drift_iv[j] = yc_iv[-1]-yc_iv[j]

# Plot trajectories
for i in range(n):
    if i >= 9:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
       
    i_sel = i
    cell_marker_i = cell_marker[i_sel]
    if i < 8:
        color_i = cmap[i]
    if i >= 8 and i < 16:
        color_i = cmap[i-8]

    #Trajectories
    # Dividing cells
    # division points; cell 3: 17, cell 5: 16, cell 7: 14, cell 9: 12
    if cell_list[i_sel][2] == 'Dividing cell':
        cell_division_point = cell_list[i_sel][3]
        
        plt.plot(cell_list[i_sel][4]+x_corr_drift_iv[-len(cell_list[i_sel][4]):], cell_list[i_sel][5]+y_corr_drift_iv[-len(cell_list[i_sel][5]):], color = 'k',
         linewidth = 3, markersize = 3) # black trajectory "daughter cell 1"
        plt.plot(cell_list[i_sel][4]+x_corr_drift_iv[-len(cell_list[i_sel][4]):], cell_list[i_sel][5]+y_corr_drift_iv[-len(cell_list[i_sel][5]):], color = color_i,
          linewidth=2.5, markersize=2.5) # colored trajectory "daughter cell 1"
        
        plt.plot(cell_list[i_sel][6]+x_corr_drift_iv[-len(cell_list[i_sel][6]):], cell_list[i_sel][7]+y_corr_drift_iv[-len(cell_list[i_sel][7]):], color = 'k',
         linewidth = 3, markersize = 3) # black trajectory "daughter cell 2"
        plt.plot(cell_list[i_sel][6]+x_corr_drift_iv[-len(cell_list[i_sel][6]):], cell_list[i_sel][7]+y_corr_drift_iv[-len(cell_list[i_sel][7]):], color = color_i,
          linewidth=2.5, markersize=2.5) # colored trajectory "daughter cell 2"
       
    # Mother cell
    # Trajectory from start to end
    plt.plot(cell_list[i_sel][0]+x_corr_drift_iv[0:len(cell_list[i_sel][0])], cell_list[i_sel][1]+y_corr_drift_iv[0:len(cell_list[i_sel][1])], color = 'k',
     linewidth = 3, markersize = 3) # line connecting start and end point black 
    plt.plot(cell_list[i_sel][0]+x_corr_drift_iv[0:len(cell_list[i_sel][0])], cell_list[i_sel][1]+y_corr_drift_iv[0:len(cell_list[i_sel][1])], color = color_i,
     linewidth=2.5, markersize=2.5) # line connecting start and end point
    
    # Dividing cells - connecting line and starting points
    if cell_list[i_sel][2] == 'Dividing cell':
        
        # Connecting line mother cell and daughter cell
        mother_x = cell_list[i_sel][0][-1]+x_corr_drift_iv[len(cell_list[i_sel][0])-1]
        mother_y = cell_list[i_sel][1][-1]+y_corr_drift_iv[len(cell_list[i_sel][1])-1]
        daughter_1_x = cell_list[i_sel][4][0]+x_corr_drift_iv[len(cell_list[i_sel][0])]
        daughter_1_y = cell_list[i_sel][5][0]+y_corr_drift_iv[len(cell_list[i_sel][1])]
        daughter_2_x = cell_list[i_sel][6][0]+x_corr_drift_iv[len(cell_list[i_sel][0])]
        daughter_2_y = cell_list[i_sel][7][0]+y_corr_drift_iv[len(cell_list[i_sel][1])]
        x_points_d1 = [mother_x, daughter_1_x]
        x_points_d2 = [mother_x, daughter_2_x]
        y_points_d1 = [mother_y, daughter_1_y]
        y_points_d2 = [mother_y, daughter_2_y]
        
        plt.plot(x_points_d1, y_points_d1, color = 'k',# linestyle = 'dotted',
         linewidth = 3, markersize = 3) # line connecting start and end point black 
        plt.plot(x_points_d1, y_points_d1, color = color_i,
         linewidth=2.5, markersize=2.5) # line connecting start and end point
        
        plt.plot(x_points_d2, y_points_d2, color = 'k', #linestyle = 'dotted',
         linewidth = 3, markersize = 3) # line connecting start and end point black 
        plt.plot(x_points_d2, y_points_d2, color = color_i,
         linewidth=2.5, markersize=2.5) # line connecting start and end point
        
        # Start point daughter cell 1
        plt.plot(cell_list[i_sel][4][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][5][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'k',
          marker='o', linewidth=20, markersize=23) # end point black edge
        plt.plot(cell_list[i_sel][4][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][5][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = color_i,
          marker='o', linewidth=20, markersize=21) # end point colored
        plt.plot(cell_list[i_sel][4][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][5][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'k',
          marker='o', linewidth=20, markersize=16) #14 end point black
        plt.plot(cell_list[i_sel][4][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][5][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'w', marker=cell_marker_i,
         linewidth=10, markersize=10*multiplying_factor) # end point number in white
        
        # Start point daughter cell 2
        plt.plot(cell_list[i_sel][6][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][7][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'k',
          marker='o', linewidth=20, markersize=23) # end point black edge
        plt.plot(cell_list[i_sel][6][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][7][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = color_i,
          marker='o', linewidth=20, markersize=21) # end point colored
        plt.plot(cell_list[i_sel][6][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][7][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'k',
          marker='o', linewidth=20, markersize=16) # end point black
        plt.plot(cell_list[i_sel][6][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][7][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'w', marker=cell_marker_i,
         linewidth=20, markersize=10*multiplying_factor) # end point number
   
    # Starting points non-dividing cells    
    plt.plot(cell_list[i_sel][0][0]+x_corr_drift_iv[0],cell_list[i_sel][1][0]+y_corr_drift_iv[0], color = 'k', marker='o',
     linewidth=2, markersize=23) # start point black edge
    plt.plot(cell_list[i_sel][0][0]+x_corr_drift_iv[0],cell_list[i_sel][1][0]+y_corr_drift_iv[0], color = color_i, marker='o',
     linewidth=2, markersize=21) # start point colored
    plt.plot(cell_list[i_sel][0][0]+x_corr_drift_iv[0],cell_list[i_sel][1][0]+y_corr_drift_iv[0], color = 'w', marker='o',
     linewidth=2, markersize=16) # start point white
    plt.plot(cell_list[i_sel][0][0]+x_corr_drift_iv[0],cell_list[i_sel][1][0]+y_corr_drift_iv[0], color = 'k', marker=cell_marker_i,
     linewidth=2, markersize=10*multiplying_factor) # start point number
    
    # End points (all cells)
    end_points = 1
    if end_points == 1:
        if cell_list[i_sel][2] == 'Dividing cell':
            # daughter cell 1
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) # end point white
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=8*multiplying_factor) # end point number
        
            # daughter cell 2
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) #13 end point white
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=10*multiplying_factor) # end point number
        
        else: # no dividing cell
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) # end point white
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=10*multiplying_factor) # end point number

#get current axes
ax = plt.gca()

#hide x-axis
ax.get_xaxis().set_visible(False)

#hide y-axis 
ax.get_yaxis().set_visible(False)

# Save figure
#plt.savefig("Mosaic_TL_mCherry/Figure_in_vitro_tracking_mosaic_icons-t31_continuous-line_black.png", dpi=600, bbox_inches='tight',pad_inches = 0)

# Plot figure
fig = plt.show()

5c. Plot (FIG 5A): correct trajectory plot for drift and rotation of gastruloid¶

In [19]:
# CORRECTION FOR DRIFT AND ROTATION
# For T = 31, plot 14 cells trajectories that are corrected for the changing rotation of the gastruloid and the 'drift' of the center's gastruloid (center of ellipse) over time.
# With drift we mean: how the center of the gastruloid - estimated by the center of the fitted ellipse - changes in position over time w.r.t. to the center of the gastruloid at the final time step.
# With rotation we mean: how the rotation of the gastruloid changes w.r.t. the gastruloid's orientation at the final time step. To correct for the rotation we use the rotated cell centers and
# rotate them again to the orientation of the gastruloid at the final time step.

%matplotlib inline
plt.rcParams['figure.figsize'] = [30, 20]

imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t31-200epochsRGB-SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])

# colors per cell
cmap = [#'#000000', #black
        '#e69d00', #orange # 1
        '#f0e442', #yellow # 2
        '#009e73', #bluish green # 3a
        '#009e73', #bluish green # 3b
        '#56b4e9', #sky blue # 4
        '#cc79a7', #reddish purple # 5a
        '#cc79a7', #reddish purple # 5b
        '#d55e00', #vermillion (reddish) # 6
        '#0072b2', #blue # 7a
        '#0072b2', #blue # 7b
        #'#ffffff', #white
        '#7f7f7f', #grey # 8
        '#e69d00', #orange # 9a
        '#e69d00', #orange # 9b
        '#f0e442', #yellow # 10
        '#009e73', #bluish green # 11
        '#56b4e9', #sky blue # 12
        '#cc79a7', #reddish purple # 13
        '#d55e00'] #vermillion (reddish) # 14

cell_marker = ['$1$','$2$','$3$','$3$','$4$','$5$','$5$','$6$','$7$','$7$','$8$','$9$','$9$','$10$','$11$','$12$','$13$','$14$']

# Rotate the rotated cell centers to the orientation of the gastruloid at the final time step
no_cells = len(rotated_cell_centers_iv)
rotated_cell_centers_rot_iv = np.empty((no_cells,t_steps,4))
for i in range(no_cells): # cells

    for j in range(t_steps): # time
        # O = position of cell
        O_x_rot_iv = np.cos(2*np.pi - theta_corr_iv[-1]) * (rotated_cell_centers_iv[i,j,0]-xc_iv[j]) + np.sin(2*np.pi - theta_corr_iv[-1]) * (rotated_cell_centers_iv[i,j,1]-yc_iv[j])
        rotated_cell_centers_rot_iv[i,j,0] = O_x_rot_iv + xc_iv[j] # rotated x-coordinate of cell
        rotated_cell_centers_rot_iv[i,j,2] = O_x_rot_iv # x component of vector from center of ellipse to O

        O_y_rot_iv = np.cos(2*np.pi - theta_corr_iv[-1]) * (rotated_cell_centers_iv[i,j,1]-yc_iv[j]) - np.sin(2*np.pi - theta_corr_iv[-1]) * (rotated_cell_centers_iv[i,j,0]-xc_iv[j])
        rotated_cell_centers_rot_iv[i,j,1] = O_y_rot_iv + yc_iv[j] # rotated y-coordinate of cell
        rotated_cell_centers_rot_iv[i,j,3] = O_y_rot_iv # y component of vector from center of ellipse to O

# Calculate drift
x_corr_drift_iv = [None] * t_steps
y_corr_drift_iv = [None] * t_steps

for j in range(0,t_steps):
    x_corr_drift_iv[j] = xc_iv[-1]-xc_iv[j]
    y_corr_drift_iv[j] = yc_iv[-1]-yc_iv[j]
    
# Plot trajectories corrected for DRIFT (add drift) + corrected for ROTATION (plot rotated-rotated cell centers)
for i in range(no_cells):
    
    # icon number size and color
    if i > 12:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
    cell_marker_i = cell_marker[i]
    color_i = cmap[i]
    
    # Plot trajectories (rotated-rotated) and add drift
    plt.plot(rotated_cell_centers_rot_iv[i,:,0]+x_corr_drift_iv, rotated_cell_centers_rot_iv[i,:,1]+y_corr_drift_iv, color = 'k', #linestyle = 'dotted',
     linewidth = 3, markersize = 3) # line connecting start and end point black 
    plt.plot(rotated_cell_centers_rot_iv[i,:,0]+x_corr_drift_iv, rotated_cell_centers_rot_iv[i,:,1]+y_corr_drift_iv, color = color_i,
     linewidth=2.5, markersize=2.5) # line connecting start and end point   
    
# Plot starting points
for i in range(no_cells):
    
    # icon number size and color
    if i > 12:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
    cell_marker_i = cell_marker[i]
    color_i = cmap[i]
    
    # For all cells, plot start icon
    plt.plot(rotated_cell_centers_rot_iv[i,0,0]+x_corr_drift_iv[0],rotated_cell_centers_rot_iv[i,0,1]+y_corr_drift_iv[0], color = 'k',
     marker='o', linewidth=20, markersize=23) # end point black edge
    plt.plot(rotated_cell_centers_rot_iv[i,0,0]+x_corr_drift_iv[0],rotated_cell_centers_rot_iv[i,0,1]+y_corr_drift_iv[0], color = color_i,
     marker='o', linewidth=20, markersize=21) # end point colored
    plt.plot(rotated_cell_centers_rot_iv[i,0,0]+x_corr_drift_iv[0],rotated_cell_centers_rot_iv[i,0,1]+y_corr_drift_iv[0], color = 'w',
     marker='o', linewidth=20, markersize=16) #14 point white
    plt.plot(rotated_cell_centers_rot_iv[i,0,0]+x_corr_drift_iv[0],rotated_cell_centers_rot_iv[i,0,1]+y_corr_drift_iv[0], color = 'k', marker=cell_marker_i,
     linewidth=10, markersize=10*multiplying_factor) # number in black
  
    # For dividing cells, only plot dividing cell (end point mother cell)
    if i == 2: # cell 3
        sel_cell = i
        div_no = cell_list[sel_cell][3] - 2# 17 > 16
    if i == 5: # cell 5
        sel_cell = i - 1
        div_no = cell_list[sel_cell][3] - 2
    if i == 8: # cell 7
        sel_cell = i - 2
        div_no = cell_list[sel_cell][3] - 2
    if i == 11: # cell 9
        sel_cell = i - 3
        div_no = cell_list[sel_cell][3] - 2

    if i == 2 or i == 5 or i == 8 or i == 11:
        #End point mother cell
        plt.plot(rotated_cell_centers_rot_iv[i,div_no,0]+x_corr_drift_iv[div_no],rotated_cell_centers_rot_iv[i,div_no,1]+y_corr_drift_iv[div_no], color = 'k',
         marker='o', linewidth=20, markersize=23) # end point black edge
        plt.plot(rotated_cell_centers_rot_iv[i,div_no,0]+x_corr_drift_iv[div_no],rotated_cell_centers_rot_iv[i,div_no,1]+y_corr_drift_iv[div_no], color = color_i,
         marker='o', linewidth=20, markersize=21) # end point colored
        plt.plot(rotated_cell_centers_rot_iv[i,div_no,0]+x_corr_drift_iv[div_no],rotated_cell_centers_rot_iv[i,div_no,1]+y_corr_drift_iv[div_no], color = 'k',
         marker='o', linewidth=20, markersize=16) #14 end point black
        plt.plot(rotated_cell_centers_rot_iv[i,div_no,0]+x_corr_drift_iv[div_no],rotated_cell_centers_rot_iv[i,div_no,1]+y_corr_drift_iv[div_no], color = 'w', marker=cell_marker_i,
         linewidth=10, markersize=10*multiplying_factor) # end point number in white
   
    # For all cells, plot end icon
    plt.plot(rotated_cell_centers_rot_iv[i,-1,0]+x_corr_drift_iv[-1],rotated_cell_centers_rot_iv[i,-1,1]+y_corr_drift_iv[-1], color = 'k',
     marker='D', linewidth=20, markersize=23) # end point black edge
    plt.plot(rotated_cell_centers_rot_iv[i,-1,0]+x_corr_drift_iv[-1],rotated_cell_centers_rot_iv[i,-1,1]+y_corr_drift_iv[-1], color = color_i,
     marker='D', linewidth=20, markersize=21) # end point colored
    plt.plot(rotated_cell_centers_rot_iv[i,-1,0]+x_corr_drift_iv[-1],rotated_cell_centers_rot_iv[i,-1,1]+y_corr_drift_iv[-1], color = 'w',
     marker='D', linewidth=20, markersize=16) #14 point white
    plt.plot(rotated_cell_centers_rot_iv[i,-1,0]+x_corr_drift_iv[-1],rotated_cell_centers_rot_iv[i,-1,1]+y_corr_drift_iv[-1], color = 'k', marker=cell_marker_i,
     linewidth=10, markersize=10*multiplying_factor) # number in black
            
#get current axes
ax = plt.gca()

#hide x-axis
ax.get_xaxis().set_visible(False)

#hide y-axis 
ax.get_yaxis().set_visible(False)

# Save figure
#plt.savefig("Mosaic_TL_mCherry/Figure_in_vitro_tracking_mosaic_t31_corrected_drift_rotation.png", dpi=600, bbox_inches='tight',pad_inches = 0)

# Plot figure
fig = plt.show()

6a. Shape analysis (FIG 5B): plot cell positions on t1 + yellow connecting lines¶

In [20]:
# T = 1, 14 cells starting points + yellow connecting lines
plt.rcParams['figure.figsize'] = [30, 20]

imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t1-200epochsRGB_SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])


cmap = [#'#000000', #black
        '#e69d00', #orange
        '#f0e442', #yellow
        '#009e73', #bluish green
        '#56b4e9', #sky blue
        '#cc79a7', #reddish purple
        '#d55e00', #vermillion (reddish)
        '#0072b2', #blue
        #'#ffffff', #white
        '#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$']
n = 14

# extra lines

#border
x_points_lines = [cell_list[6][0][0],cell_list[2][0][0],cell_list[1][0][0],cell_list[0][0][0],
                  cell_list[5][0][0],cell_list[11][0][0],cell_list[13][0][0],cell_list[12][0][0],cell_list[9][0][0],cell_list[6][4][0]]
y_points_lines = [cell_list[6][1][0],cell_list[2][1][0],cell_list[1][1][0],cell_list[0][1][0],
                  cell_list[5][1][0],cell_list[11][1][0],cell_list[13][1][0],cell_list[12][1][0],cell_list[9][1][0],cell_list[6][1][0]]
#in between
x_points_lines_mid_1 = [cell_list[2][0][0],cell_list[3][0][0],cell_list[5][0][0]]
y_points_lines_mid_1 = [cell_list[2][1][0],cell_list[3][1][0],cell_list[5][1][0]]
x_points_lines_mid_2 = [cell_list[6][0][0],cell_list[7][0][0],cell_list[5][0][0]]
y_points_lines_mid_2 = [cell_list[6][1][0],cell_list[7][1][0],cell_list[5][1][0]]
x_points_lines_mid_3 = [cell_list[9][0][0],cell_list[10][0][0],cell_list[11][0][0]]
y_points_lines_mid_3 = [cell_list[9][1][0],cell_list[10][1][0],cell_list[11][1][0]]

plt.plot(x_points_lines, y_points_lines,'y',linewidth = 5)
plt.plot(x_points_lines_mid_1, y_points_lines_mid_1,'y',linewidth = 5)
plt.plot(x_points_lines_mid_2, y_points_lines_mid_2,'y',linewidth = 5)
plt.plot(x_points_lines_mid_3, y_points_lines_mid_3,'y',linewidth = 5)

for i in range(n):
    if i >= 9:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
       
    i_sel = i
    cell_marker_i = cell_marker[i_sel]
    if i < 8:
        color_i = cmap[i]
    if i >= 8 and i < 16:
        color_i = cmap[i-8]

    # start points   
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker='o',
     linewidth=2, markersize=23) # start point black edge
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = color_i, marker='o',
     linewidth=2, markersize=21) # start point colored
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'w', marker='o',
     linewidth=2, markersize=16) # start point white
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker=cell_marker_i,
     linewidth=2, markersize=10*multiplying_factor) # start point number
     
#get current axes
ax = plt.gca()

#hide x-axis
ax.get_xaxis().set_visible(False)

#hide y-axis 
ax.get_yaxis().set_visible(False)

# Save figure
#plt.savefig("Mosaic_TL_mCherry/Figure_in_vitro_tracking_mosaic_icons_lines_t1_more_white.png", dpi=600, bbox_inches='tight',pad_inches = 0)

# Plot figure
fig = plt.show()

6b. Shape analysis (FIG 5B): plot cell positions on t31 + yellow connecting lines¶

In [21]:
# T = 31, 14 cells trajectories
%matplotlib inline
plt.rcParams['figure.figsize'] = [30, 20]

imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t31-200epochsRGB-SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])


cmap = [#'#000000', #black
        '#e69d00', #orange
        '#f0e442', #yellow
        '#009e73', #bluish green
        '#56b4e9', #sky blue
        '#cc79a7', #reddish purple
        '#d55e00', #vermillion (reddish)
        '#0072b2', #blue
        #'#ffffff', #white
        '#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$']
n = 14

# extra lines

#border
x_points_lines = [cell_list[6][4][-1],cell_list[2][4][-1],cell_list[1][0][-1],cell_list[0][0][-1],
                  cell_list[5][0][-1],cell_list[11][0][-1],cell_list[13][0][-1],cell_list[12][0][-1],cell_list[9][0][-1],cell_list[6][4][-1]]
y_points_lines = [cell_list[6][5][-1],cell_list[2][5][-1],cell_list[1][1][-1],cell_list[0][1][-1],
                  cell_list[5][1][-1],cell_list[11][1][-1],cell_list[13][1][-1],cell_list[12][1][-1],cell_list[9][1][-1],cell_list[6][5][-1]]
#in between
x_points_lines_mid_1 = [cell_list[2][4][-1],cell_list[3][0][-1],cell_list[5][0][-1]]
y_points_lines_mid_1 = [cell_list[2][5][-1],cell_list[3][1][-1],cell_list[5][1][-1]]
x_points_lines_mid_2 = [cell_list[6][4][-1],cell_list[7][0][-1],cell_list[5][0][-1]]
y_points_lines_mid_2 = [cell_list[6][5][-1],cell_list[7][1][-1],cell_list[5][1][-1]]
x_points_lines_mid_3 = [cell_list[9][0][-1],cell_list[10][0][-1],cell_list[11][0][-1]]
y_points_lines_mid_3 = [cell_list[9][1][-1],cell_list[10][1][-1],cell_list[11][1][-1]]

plt.plot(x_points_lines, y_points_lines,'y',linewidth = 5)
plt.plot(x_points_lines_mid_1, y_points_lines_mid_1,'y',linewidth = 5)
plt.plot(x_points_lines_mid_2, y_points_lines_mid_2,'y',linewidth = 5)
plt.plot(x_points_lines_mid_3, y_points_lines_mid_3,'y',linewidth = 5)

for i in range(n):
    if i >= 9:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
       
    i_sel = i
    cell_marker_i = cell_marker[i_sel]
    if i < 8:
        color_i = cmap[i]
    if i >= 8 and i < 16:
        color_i = cmap[i-8]

    # end points
    end_points = 1
    if end_points == 1:
        if cell_list[i_sel][2] == 'Dividing cell':
            # daughter cell 1
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) # end point white
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=8*multiplying_factor) # end point number
        
            # daughter cell 2
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) # end point white
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=10*multiplying_factor) # end point number
        
        else: # no dividing cell
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) # end point white
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=10*multiplying_factor) # end point number
            
#get current axes
ax = plt.gca()

#hide x-axis
ax.get_xaxis().set_visible(False)

#hide y-axis 
ax.get_yaxis().set_visible(False)

# Save figure
#plt.savefig("Figure_in_vitro_tracking_mosaic_icons_lines_t31_endpointsonly.png", dpi=600, bbox_inches='tight',pad_inches = 0)

# Plot figure
fig = plt.show()

7a. Analysis/plot: calculate the distance of cells moved along the long axis and plot for final position on long axis¶

In [22]:
t_end_iv = 30 # t = 31
t_start_iv = 0 # t = 1

# Determine the length of the gastruloid at t_end_iv = 2 * long arm of ellipse fit. The long arm should be stored in the
# parameter 'b', and the short arm in parameter 'a', but the fitting tool sometimes stored them incorrectly.
# Therefore, check the value of a and b first to determine which parameter corresponds to the long arm.
# In our in vitro example, the long arm of the fitted ellipses should be on the y-axis. If a and b are incorrectly switched,
# the values stored in the x coordinates correspond to the long arm.

if a_iv[t_end_iv] > b_iv[t_end_iv]:
    
    length_gastruloid_iv = 2*a_iv[t_end_iv] # used to normalize the distance the cells moved along the long axis
    print('a > b, so we use the rotated x-coordinates to calculate the distance cells moved along the long axis of the gastruloid')
     
    end_pos_on_long_axis_iv = rotated_cell_centers_iv[:,t_end_iv,0] # rotated x-coordinate at t = t_end_iv
    end_pos_on_long_axis_iv_norm = (end_pos_on_long_axis_iv - xc_iv[t_end_iv])/(length_gastruloid_iv) # (rotated x-coordinate at t_end_iv - x_center_of_ellipse at t_end_iv) divided by length gastruloid
    moved_end_iv_norm = ((end_pos_on_long_axis_iv - rotated_cell_centers_iv[:,t_start_iv,0]) - (xc_iv[t_end_iv] - xc_iv[t_start_iv]) ) / length_gastruloid_iv # to normalize, divide by length of gastruloid

else: 
    length_gastruloid_iv = 2*b_iv[t_end_iv] # used to normalize the distance the cells moved along the long axis
    print('a < b, so we use the rotated y-coordinates to calculate the distance cells moved along the long axis of the gastruloid')
      
    end_pos_on_long_axis_iv = rotated_cell_centers_iv[:,t_end_iv,1] # rotated y-coordinate at t = t_end_iv
    end_pos_on_long_axis_iv_norm = (end_pos_on_long_axis_iv - yc_iv[t_end_iv])/(length_gastruloid_iv) # (rotated y-coordinate at t_end_iv - y_center_of_ellipse at t_end_iv) divided by length gastruloid
    moved_end_iv_norm = ((end_pos_on_long_axis_iv - rotated_cell_centers_iv[:,t_start_iv,1]) - (yc_iv[t_end_iv] - yc_iv[t_start_iv]) ) / length_gastruloid_iv # to normalize, divide by length of gastruloid
    
# -------- FIGURE ----------
# Figure parameters    
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'  
plt.xlabel('End position projected on long axis', fontsize = 16)
plt.ylabel('Movement along long axis', fontsize = 16)
plt.xlim(-0.55, 0.55)
plt.ylim(-0.51, 0.51)
plt.xticks([-0.5, -0.25, 0, 0.25, 0.5],['-L/2','-L/4','0','L/4','L/2'])
plt.yticks([-0.5, -0.25, 0, 0.25, 0.5],['-L/2','-L/4','0','L/4','L/2'])

# In vitro time lapse data
ax.scatter(end_pos_on_long_axis_iv_norm, moved_end_iv_norm, color = '#cc79a7', marker = 'o', label = 'in vitro',s = 50)

handles, labels = ax.get_legend_handles_labels()
lgd = ax.legend(handles, labels, loc='upper left', bbox_to_anchor=(1,1))
fig.set_size_inches((110/5)/2.54, (76/5)/2.54)

plt.show()
a < b, so we use the rotated y-coordinates to calculate the distance cells moved along the long axis of the gastruloid

7b. Analysis/plot: calculate the distance moved along short axis and plot for final position on long axis¶

In [23]:
t_end_iv = 30 # t = 31
t_start_iv = 0 # t = 1

# Determine the length of the gastruloid at t_end_iv = 2 * long arm of ellipse fit. The long arm should be stored in the
# parameter 'b', and the short arm in parameter 'a', but the fitting tool sometimes stored them incorrectly.
# Therefore, check the value of a and b first to determine which parameter corresponds to the long or short arm.
# In our in vitro example, the long arm of the fitted ellipses should be on the y-axis. The short arm on the x-axis. 
# If a and b are incorrectly switched, the values stored in the x coordinates correspond to the long arm, and y coordinates correspond to the short arm.

if a_iv[t_end_iv] > b_iv[t_end_iv]:
    
    length_gastruloid_iv_short = 2*b_iv[t_end_iv] # used to normalize the distance the cells moved along the short axis
    length_gastruloid_iv_long = 2*a_iv[t_end_iv] # used to normalize the position of the cells along the long axis
    print('a > b, so we use the rotated y-coordinates to calculate the distance cells moved along the short axis of the gastruloid')
      
    end_pos_on_short_axis_iv = rotated_cell_centers_iv[:,t_end_iv,1] # rotated y-coordinate at t = t_end_iv
    end_pos_on_long_axis_iv = rotated_cell_centers_iv[:,t_end_iv,0] # rotated x-coordinate at t = t_end_iv
    end_pos_on_short_axis_iv_norm = (end_pos_on_short_axis_iv - yc_iv[t_end_iv])/(length_gastruloid_iv_short) # (rotated y-coordinate at t_end_iv - y_center_of_ellipse at t_end_iv) divided by length gastruloid
    end_pos_on_long_axis_iv_norm = (end_pos_on_long_axis_iv - xc_iv[t_end_iv])/(length_gastruloid_iv_long) # (rotated y-coordinate at t_end_iv - y_center_of_ellipse at t_end_iv) divided by length gastruloid

    moved_end_iv_short_norm = ((end_pos_on_short_axis_iv - rotated_cell_centers_iv[:,t_start_iv,1]) - (yc_iv[t_end_iv] - yc_iv[t_start_iv]) ) / length_gastruloid_iv_short # to normalize, divide by length of gastruloid
    moved_end_iv_long_norm = ((end_pos_on_long_axis_iv - rotated_cell_centers_iv[:,t_start_iv,0]) - (xc_iv[t_end_iv] - xc_iv[t_start_iv]) ) / length_gastruloid_iv_long # to normalize, divide by length of gastruloid

else: 
    length_gastruloid_iv_short = 2*a_iv[t_end_iv] # used to normalize the distance the cells moved along the short axis
    length_gastruloid_iv_long = 2*b_iv[t_end_iv] # # used to normalize the position of the cells along the long axis
    print('a < b, so we use the rotated x-coordinates to calculate the distance cells moved along the short axis of the gastruloid')
     
    end_pos_on_short_axis_iv = rotated_cell_centers_iv[:,t_end_iv,0] # rotated x-coordinate at t = t_end_iv
    end_pos_on_long_axis_iv = rotated_cell_centers_iv[:,t_end_iv,1] # rotated y-coordinate at t = t_end_iv
    end_pos_on_short_axis_iv_norm = (end_pos_on_short_axis_iv - xc_iv[t_end_iv])/(length_gastruloid_iv_short) # (rotated x-coordinate at t_end_iv - x_center_of_ellipse at t_end_iv) divided by length gastruloid
    end_pos_on_long_axis_iv_norm = (end_pos_on_long_axis_iv - yc_iv[t_end_iv])/(length_gastruloid_iv_long) # (rotated y-coordinate at t_end_iv - y_center_of_ellipse at t_end_iv) divided by length gastruloid

    moved_end_iv_short_norm = ((end_pos_on_short_axis_iv - rotated_cell_centers_iv[:,t_start_iv,0]) - (xc_iv[t_end_iv] - xc_iv[t_start_iv]) ) / length_gastruloid_iv_short # to normalize, divide by length of gastruloid
    moved_end_iv_long_norm = ((end_pos_on_long_axis_iv - rotated_cell_centers_iv[:,t_start_iv,1]) - (yc_iv[t_end_iv] - yc_iv[t_start_iv]) ) / length_gastruloid_iv_long # to normalize, divide by length of gastruloid

# -------- FIGURE ----------
# Figure parameters    
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'  
plt.xlabel('End position projected on short axis', fontsize = 16)
plt.ylabel('Movement along short axis', fontsize = 16)
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.xticks([-1, -0.5, 0, 0.5, 1],['S','-S/2','0','S/2', 'S'])
plt.yticks([-1,-0.5, 0, 0.5, 1],['-S','-S/2','0','S/2','S'])

ax.scatter(end_pos_on_short_axis_iv_norm, moved_end_iv_short_norm, color = '#cc79a7', marker = 'o', s = 50)#label = 'in vitro',s = 50)

fig.set_size_inches((110/5)/2.54, (76/5)/2.54)

plt.show()
a < b, so we use the rotated x-coordinates to calculate the distance cells moved along the short axis of the gastruloid

In silico¶

Create a binary image¶

In [24]:
# Read image
%matplotlib inline
img = cv2.imread('Mosaic_3/image0099900.png',0)
plt.imshow(img,cmap='gray')
plt.show()
In [25]:
# Erode the image multiple times, and subsequently dilate the image multiple times

%matplotlib inline
# Taking a matrix of size 5 as the kernel
kernel = np.ones((5, 5), np.uint8)

# The first parameter is the original image,
# kernel is the matrix with which image is
# convolved and third parameter is the number
# of iterations, which will determine how much
# you want to erode/dilate a given image.

img_erosion = cv2.erode(img, kernel, iterations=4)
img_dilation = cv2.dilate(img_erosion, kernel, iterations=4)

imagem = cv2.bitwise_not(img_erosion)

plt.imshow(imagem,cmap='gray')
plt.show()
plt.imshow(img_dilation,cmap='gray')
plt.show()
# save image
#cv2.imwrite('Mosaic_3/image0099900-binary.png', imagem)

Import images and create a GIF¶

In [ ]:
# read the images and store them in a list in the correct order
image_path = Path('Mosaic_corrected')
#images = list(image_path.glob('*.png'))
image_list = []
for i in range(0,100000,100):
    n = i
    STRT = "Mosaic_corrected/image"
    EXTN = ".png"
    filename = STRT + str(n).zfill(7) + EXTN
    image_list.append(imageio.imread(filename))
In [ ]:
# from the list of images, write a GIF (standard 10 fps)
imageio.mimwrite('Mosaic_corrected/animated_mosaic_corrected.gif',image_list)

Alternatively, one coud also open the images in ImageJ, convert it to 8-bit, add a Monte Carlo step time stamp and save the sequence as a GIF.

In silico tracking - mosaic example¶

0a. Create mask and extract outline¶

In [27]:
def extract_outline():

    for x in range(0,1000):
        ind = list()
        cont_all = list()
        image_loc = 'Mosaic_corrected/image'+"{:05d}".format(x)+'00.png'
        file_loc = 'Mosaic_corrected/outlines/outline_'+str(x) + '.txt'
        #file1 = open(file_loc,"a")
        mask = cv2.imread(image_loc)
        mask = 255-mask
        mask = mask[1:-1,1:-1]
        mask = mask[1:-1,1
                    :-1]
        mask = mask[1:-1,1:-1]

        if mask is not None:

            imgray = cv2.cvtColor(mask,cv2.COLOR_BGR2GRAY)
            ret,thresh = cv2.threshold(imgray,200,255,0)
            contours, im2 = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_NONE)         
            maxsizeloc = 0
            maxsize = 0
            for i in range (len(contours)):
                if contours[i].shape[0] > maxsize:
                    maxsizeloc = i
                    maxsize = contours[i].shape[0]
            cont = np.zeros((contours[maxsizeloc].shape[0], 2))
            cont[:,0] = contours[maxsizeloc][:,0,0]
            cont[:,1] = contours[maxsizeloc][:,0,1]

            plt.plot(cont[:,0], cont[:,1])
            cont_all.append(cont)
            
            #for i in range (maxsize):
                #file1.write(str(cont_all[0][i][0]) + ',' + str(cont_all[0][i][1]) + '\n')
            #file1.close()

    return(cont_all)

extract_outline()
Out[27]:
[array([[420., 296.],
        [420., 297.],
        [419., 298.],
        ...,
        [423., 296.],
        [422., 296.],
        [421., 296.]])]

0b. Load outline and fit ellipse, store ellipse parameters¶

In [28]:
# Ellipse fit parameters
xc = np.empty(1000) # center coordinate x of ellipse
yc = np.empty(1000) # center coordinate y of ellipse
a = np.empty(1000) # length short arm of ellipse
b = np.empty(1000) # length long arm of ellipse
theta= np.empty(1000) # rotation angle of ellipse

# Figure parameters
%matplotlib inline
fig, axs = plt.subplots(figsize = (10,10))#, sharex=True, sharey=True)#, figsize=(5,15))
cmap = [#'#000000', #black
        '#e69d00', #orange
        '#56b4e9', #sky blue
        '#009e73', #bluish green
        '#f0e442', #yellow
        '#0072b2', #blue
        '#d55e00', #vermillion (reddish)
        '#cc79a7', #reddish purple
        #'#ffffff', #white
        '#7f7f7f'] #grey

for i in range(1000):
    # Load outline coordinates for MCS i
    with open("Mosaic_corrected/outline_" + str(i) + ".txt") as file_name:
        points = np.loadtxt(file_name, delimiter=",")
    a_points = np.array(points)
    x = a_points[:, 0]
    y = a_points[:, 1]
    
    #Calculate ellipse fit from outline and store parameters
    ell = EllipseModel()
    ell.estimate(a_points)
    xc[i], yc[i], a[i], b[i], theta[i] = ell.params
    
    # Write parameters to txt file
    #file_loc = 'Mosaic_corrected/ellipse_params/MCS_'+str(i)+'-x_y_theta_a_b.txt'
    #file1 = open(file_loc,"a")
    #file1.write(str(xc[i]) + ',' + str(yc[i]) + ',' + str(theta[i]) + ',' + str(a[i]) + ',' + str(b[i]) + '\n')
    
    # Plot ellipse fit for one MCS
    t_end = 999
    if i == t_end:
        axs.scatter(x, y, color = cmap[4])
        axs.set_ylim(0,1000)
        axs.set_xlim(0,1000)
        axs.scatter(xc[i], yc[i], color=cmap[0], s=50)
        ell_patch = Ellipse((xc[i], yc[i]), 2*a[i], 2*b[i], theta[i]*180/np.pi, edgecolor=cmap[0], facecolor='none', linewidth = 3)
        axs.add_patch(ell_patch)

        if a[i] < b[i]:
            print('a',a[i],' < b',b[i], 'theta = ', theta[i]*180/np.pi)

# Save figure
#plt.savefig("Mosaic_corrected/ellipse_fit_t999.png", dpi=600, pad_inches = 0)    
    
plt.show()
print('Theta at MCS 100000 = ', theta[999]*180/np.pi)
Theta at MCS 100000 =  53.37818613761889
In [29]:
# Plot length of long arm and short arm over MCS - use to determine range of MCS to plot
%matplotlib inline
x1000 = list(range(0,1000))
for i in x1000:
    if a[i] > b[i]: # a is long arm, b is short arm
        plt.scatter(x1000[i],a[i], s = 0.5, color = 'k')
        plt.scatter(x1000[i],b[i], s = 0.5, color = 'b')
    elif a[i] < b[i]: # b is long arm, a is short arm
        plt.scatter(x1000[i],b[i], s = 0.5, color = 'k')
        plt.scatter(x1000[i],a[i], s = 0.5, color = 'b') 
t_start = 199
line_1 = t_start
plt.plot([line_1, line_1],[0, 220])
plt.show()

1a. Store x- and y-coordinates of selected cells (20) for selected MCS (16) in arrays¶

In [30]:
# Define number of cells and the selected time points you want to plot the trajectory for
No_cells = 20 # number of cells that are being tracked
selected_points = [t_start,249,299,349,399,449,499,549,599,649,699,749,799,849,899,949,999] # 16 steps
    
# Define list of arrays for x- and y-coordinates for each tracked cell
path_x_sim = [None] * (No_cells)
path_y_sim = [None] * (No_cells)
path_x_sim_final = [None] * No_cells
path_y_sim_final = [None] * No_cells
In [31]:
# Store x- and y-coordinates for each tracked cell in lists 'path_x_sim_final' and 'path_y_sim_final'
for i in range(10,int(No_cells*10+10),10):
    filepath = ""
    file_dir = ['Mosaic_corrected/cell_','.txt']
    filepath = file_dir[0]+str(i)+file_dir[1]
    with open(filepath) as f:
            for line in f:
                currentline = line.split(",")
                x_line = float(currentline[0])
                y_line = float(currentline[1])
                ii = int(i/10-1)
                path_x_sim[ii] = np.append(path_x_sim[ii],x_line)
                path_y_sim[ii] = np.append(path_y_sim[ii],y_line)
    path_x_sim_final[ii] = path_x_sim[ii][1:(len(path_x_sim[ii])+1)]
    path_y_sim_final[ii] = path_y_sim[ii][1:(len(path_x_sim[ii])+1)]
In [32]:
# Store x- and y-coordinates of each tracked cell for only the selected time points in lists 'x_sel_final' and 'y_sel_final'
x_sel = [None] * No_cells
y_sel = [None] * No_cells
x_sel_final = [None] * No_cells
y_sel_final = [None] * No_cells

for j in range(No_cells):
    for k in range(len(selected_points)):
        x_sel[j] = np.append(x_sel[j],2*path_x_sim_final[j][selected_points[k]])
        y_sel[j] = np.append(y_sel[j],2*path_y_sim_final[j][selected_points[k]])
    x_sel_final[j] = x_sel[j][1:(len(selected_points)+1)]
    y_sel_final[j] = y_sel[j][1:(len(selected_points)+1)]

2. Plot: 12 start positions on start image¶

In [33]:
%matplotlib inline

# Select image (time point) and adjust size (optional)
imgplot_100= plt.imread('Mosaic_corrected/image0019900.png')
plt.rcParams['figure.figsize'] = [200, 100]

# Make image gray (so no dark blue cells)
R_100, G_100, B_100 = imgplot_100[:,:,0], imgplot_100[:,:,1], imgplot_100[:,:,2]
imgGray_100 = 0.1 * R_100 + 0.1 * G_100 + 0.7*B_100

cmap = [#'#000000', #black
        '#e69d00', #orange
        '#56b4e9', #sky blue
        '#009e73', #bluish green
        '#f0e442', #yellow
        '#0072b2', #blue
        '#d55e00', #vermillion (reddish)
        '#7f7f7f', #grey
        '#cc79a7']#, #reddish purple
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$','$15$','$16$','$17$','$18$','$19$','$20$']
markers_numbered_selected = [11,10,14,17,16,18,3,5,6,7,8,9] # 12 cells

for i in range(len(markers_numbered_selected)):
    
    if i <= 5:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
    
    i_sel_rot = (markers_numbered_selected[i])*10-1
    i_sel = markers_numbered_selected[i]-1
    cell_marker_i = cell_marker[i_sel]
    if i == 0:
        color_i = cmap[2]
    if i == 2:
        color_i = cmap[7]
    if i > 0 and i < 8 and i != 2:
        color_i = cmap[i]
    if i == 8:
        color_i = cmap[0]
    if i == 9:
        color_i = cmap[2]
    if i == 10:
        color_i = cmap[1]
    if i == 11:
        color_i = cmap[3] 
    
    t_start_sel = 0
    plt.plot(x_sel_final[i_sel][t_start_sel],y_sel_final[i_sel][t_start_sel], color = 'k', marker='o',
     linewidth=2, markersize=130) # start point black edge
    plt.plot(x_sel_final[i_sel][t_start_sel],y_sel_final[i_sel][t_start_sel], color = color_i, marker='o',
     linewidth=2, markersize=120) # start point colored
    plt.plot(x_sel_final[i_sel][t_start_sel],y_sel_final[i_sel][t_start_sel], color = 'w', marker='o',
     linewidth=2, markersize=80) # start point white
    plt.plot(x_sel_final[i_sel][t_start_sel],y_sel_final[i_sel][t_start_sel], color = 'k', marker=cell_marker_i,
     linewidth=20, markersize=40*multiplying_factor) # start point number

# Save figure
plt.imshow(imgGray_100, cmap='gray')
# Save figure
#plt.savefig("Mosaic_corrected/Figure_in_silico_mosaic_corrected_12_icons_t199.png", bbox_inches='tight',pad_inches = 0)
plt.show()

3. Plot: and 12 end positions on end image¶

In [34]:
%matplotlib inline

# Select image (time point) and adjust size (optional)
imgplot= plt.imread('Mosaic_corrected/image0099900.png') # final time point
plt.rcParams['figure.figsize'] = [200, 100]

# Make image gray (so no dark blue cells)
R, G, B = imgplot[:,:,0], imgplot[:,:,1], imgplot[:,:,2]
imgGray = 0.1 * R + 0.1 * G + 0.7*B #0.2989 * R + 0.5870 * G + 0.1140 * B

cmap = [#'#000000', #black
        '#e69d00', #orange
        '#56b4e9', #sky blue
        '#009e73', #bluish green
        '#f0e442', #yellow
        '#0072b2', #blue
        '#d55e00', #vermillion (reddish)
        '#7f7f7f', #grey
        '#cc79a7']#, #reddish purple
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$','$15$','$16$','$17$','$18$','$19$','$20$']
markers_numbered_selected = [11,10,14,17,16,18,3,5,6,7,8,9]

for i in range(len(markers_numbered_selected)):
    
    if i <= 5:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
    
    i_sel_rot = (markers_numbered_selected[i])*10-1
    i_sel = markers_numbered_selected[i]-1
    cell_marker_i = cell_marker[i_sel]
    if i == 0:
        color_i = cmap[2]
    if i == 2:
        color_i = cmap[7]
    if i > 0 and i < 8 and i != 2:
        color_i = cmap[i]
    if i == 8:
        color_i = cmap[0]
    if i == 9:
        color_i = cmap[2]
    if i == 10:
        color_i = cmap[1]
    if i == 11:
        color_i = cmap[3] 
   
    plt.plot(x_sel_final[i_sel][-1],y_sel_final[i_sel][-1], color = 'k',
      marker='D', linewidth=20, markersize=130) # end point black edge
    plt.plot(x_sel_final[i_sel][-1],y_sel_final[i_sel][-1], color = color_i,
      marker='D', linewidth=20, markersize=120) # end point colored
    plt.plot(x_sel_final[i_sel][-1],y_sel_final[i_sel][-1], color = 'w',
      marker='D', linewidth=20, markersize=70) # end point white
    plt.plot(x_sel_final[i_sel][-1],y_sel_final[i_sel][-1], color = 'k', marker=cell_marker_i,
     linewidth=20, markersize=40*multiplying_factor) # end point number
    
plt.imshow(imgGray, cmap='gray')  
plt.show()

4. [Trajectories preparation] Rotate cell centers with found theta from ellipse fit, rotate around ellipse center¶

In [36]:
# Read in cell centers coordinates and store in cell_centers[cell][time][coordinate]
cell_centers = []
for i in range(1,201):
    with open("Mosaic_corrected/cell_" + str(i) + ".txt") as file_name:
        cell_centers.append(np.loadtxt(file_name, delimiter=","))

for cell in range(200):
    for time in range(1000):
        for coordinate in range(2):
            cell_centers[cell][time][coordinate] = 2*cell_centers[cell][time][coordinate]

# Rotate cell centers and store in rotated_cell_centers (rotated around ellipse center)
rotated_cell_centers = np.empty((200, 1000, 4))
theta_corr = np.empty(len(theta))

for i in range(1,201): # cells
  
    for j in range(len(cell_centers[0])): # time / MCS
         
        # correct theta if it shows a jump (= orientation error)  
        if theta[j]*180/np.pi >64 : # gastruloid angle
            theta_corr[j] = theta[j] - 1.57079633
        elif theta[j]*180/np.pi <0 : 
            theta_corr[j] = theta[j] + 1.57079633
        else:
            theta_corr[j] = theta[j]

        # Calculate rotated coordinates, rotate around center of ellips
        # xc/yc = center coordinates of ellipse, O_x/O_y = coordinates of cell
        O_x = np.cos(theta_corr[j]) * (cell_centers[i-1][j][0]-xc[j]) + np.sin(theta_corr[j]) * (cell_centers[i-1][j][1]-yc[j])
        rotated_cell_centers[i-1,j,0] = O_x + xc[j] # rotated x-coordinate of cell
        rotated_cell_centers[i-1,j,2] = O_x # x component of vector from center of ellipse to position of cell
        O_y = np.cos(theta_corr[j]) * (cell_centers[i-1][j][1]-yc[j]) - np.sin(theta_corr[j]) * (cell_centers[i-1][j][0]-xc[j])
        rotated_cell_centers[i-1,j,1] = O_y + yc[j] # rotated y-coordinate of cell
        rotated_cell_centers[i-1,j,3] = O_y # y component of vector from center of ellipse to position of cell  
In [37]:
# Investigate theta_corr
%matplotlib inline
print('Theta at t_start: ', theta[t_start]*180/np.pi, 'and theta at t_end: ',theta[999]*180/np.pi)
print('Corrected theta at t_start: ', theta_corr[t_start]*180/np.pi, 'and corrected theta at t_end: ',theta_corr[999]*180/np.pi)
theta_x = range(0,1000,1)
# Plot theta
plt.scatter(theta_x,theta*180/np.pi, s = 1, color = 'red')
# Plot theta_corr
plt.scatter(theta_x,theta_corr*180/np.pi, s = 1, color = 'black')

# Plot t_start
plt.scatter(199,theta_corr[199]*180/np.pi, s = 30, color = 'limegreen')
plt.xlim(-10, 1010)
Theta at t_start:  127.77361178960888 and theta at t_end:  53.37818613761889
Corrected theta at t_start:  37.773611605969975 and corrected theta at t_end:  53.37818613761889
Out[37]:
(-10.0, 1010.0)

5. Plot: trajectories on end image (+ correction for drift + correction for rotation)¶

In [38]:
# CORRECTION FOR DRIFT AND FOR ROTATION
# Rotate rotated cell centers back to orientation gastruloid at final MCS, plot trajectories at final MCS and correct for drift
# Plot figure: each trajectory a different color (max. 8 colors), color in the color-blind-friendly palette

%matplotlib inline

# Select image (time point) and adjust size (optional)
imgplot= plt.imread('Mosaic_corrected/image0099900.png') # final time point
plt.rcParams['figure.figsize'] = [200, 100]
# Make image gray (so no dark blue cells)
R, G, B = imgplot[:,:,0], imgplot[:,:,1], imgplot[:,:,2]
imgGray = 0.1 * R + 0.1 * G + 0.7*B #0.2989 * R + 0.5870 * G + 0.1140 * B

# Calculate drift of gastruloid for each MCS w.r.t. final MCS using the difference in position of ellipse centers
x_corr_drift_all = np.empty((1000))
y_corr_drift_all = np.empty((1000))
# Rotate all rotated cell centers back to gastruloid's orientation at the final MCS (999)
rotated_cell_centers_rot = np.empty((200, 1000, 4))

cmap = [#'#000000', #black
        '#e69d00', #orange
        '#56b4e9', #sky blue
        '#009e73', #bluish green
        '#f0e442', #yellow
        '#0072b2', #blue
        '#d55e00', #vermillion (reddish)
        '#7f7f7f', #grey
        '#cc79a7']#, #reddish purple
        #'#ffffff', #white
        #'#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$','$15$','$16$','$17$','$18$','$19$','$20$']
markers_numbered_selected = [11,10,14,17,16,18,3,5,6,7,8,9]

for i in range(0,200): # cells
    
    for j in range(0,1000): # MCSs
        x_corr_drift_all[j] = xc[-1]-xc[j]
        y_corr_drift_all[j] = yc[-1]-yc[j]
        
        O_x_rot = np.cos(2*np.pi - theta_corr[-1]) * (rotated_cell_centers[i][j][0]-xc[j]) + np.sin(2*np.pi - theta_corr[-1]) * (rotated_cell_centers[i][j][1]-yc[j])
        rotated_cell_centers_rot[i,j,0] = O_x_rot + xc[j] # rotated x-coordinate of cell
        rotated_cell_centers_rot[i,j,2] = O_x_rot # x component of vector from center of ellipse to position of cell
        
        O_y_rot = np.cos(2*np.pi - theta_corr[-1]) * (rotated_cell_centers[i][j][1]-yc[j]) - np.sin(2*np.pi - theta_corr[-1]) * (rotated_cell_centers[i][j][0]-xc[j])
        rotated_cell_centers_rot[i,j,1] = O_y_rot + yc[j] # rotated y-coordinate of cell
        rotated_cell_centers_rot[i,j,3] = O_y_rot # y component of vector from center of ellipse to position of cell

for i in range(len(markers_numbered_selected)): # loop over cells
    i_sel_rot = (markers_numbered_selected[i])*10-1
    i_sel = (markers_numbered_selected[i])-1
    cell_marker_i = cell_marker[i_sel]
    
    if i == 0:
        color_i = cmap[2]
    if i == 2:
        color_i = cmap[7]
    if i > 0 and i < 8 and i != 2:
        color_i = cmap[i]
    if i == 8:
        color_i = cmap[0]
    if i == 9:
        color_i = cmap[2]
    if i == 10:
        color_i = cmap[1]
    if i == 11:
        color_i = cmap[3]
    
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points,0] + x_corr_drift_all[selected_points], rotated_cell_centers_rot[i_sel_rot,selected_points,1] + y_corr_drift_all[selected_points], color = 'k',
             linewidth=25, markersize=25) # line connecting start and end point black
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points,0] + x_corr_drift_all[selected_points], rotated_cell_centers_rot[i_sel_rot,selected_points,1] + y_corr_drift_all[selected_points], color = color_i,
             linewidth=20, markersize=20) # line connecting start and end point, linewidth=25, markersize=25

for i in range(len(markers_numbered_selected)):
    
    if i <= 5:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
    
    i_sel_rot = (markers_numbered_selected[i])*10-1
    i_sel = markers_numbered_selected[i]-1 #i
    cell_marker_i = cell_marker[i_sel]
    if i == 0:
        color_i = cmap[2]
    if i == 2:
        color_i = cmap[7]
    if i > 0 and i < 8 and i != 2:
        color_i = cmap[i]
    if i == 8:
        color_i = cmap[0]
    if i == 9:
        color_i = cmap[2]
    if i == 10:
        color_i = cmap[1]
    if i == 11:
        color_i = cmap[3]
    
    # Calculate drift for starting point w.r.t. the final point and add that to starting positions
    x_corr_drift_f = xc[999]-xc[t_start]
    y_corr_drift_f = yc[999]-yc[t_start]
    # plot start point
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[0],0] + x_corr_drift_all[selected_points[0]], rotated_cell_centers_rot[i_sel_rot,selected_points[0],1] + y_corr_drift_all[selected_points[0]],
     color = 'k', marker='o',linewidth=2, markersize=130) # start point black edge
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[0],0] + x_corr_drift_all[selected_points[0]], rotated_cell_centers_rot[i_sel_rot,selected_points[0],1] + y_corr_drift_all[selected_points[0]],
     color = color_i, marker='o', linewidth=2, markersize=120) # start point colored,  use alpha = .5 for transparency
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[0],0] + x_corr_drift_all[selected_points[0]], rotated_cell_centers_rot[i_sel_rot,selected_points[0],1] + y_corr_drift_all[selected_points[0]],
     color = 'w', marker='o', linewidth=2, markersize=80)#70) # start point white
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[0],0] + x_corr_drift_all[selected_points[0]], rotated_cell_centers_rot[i_sel_rot,selected_points[0],1] + y_corr_drift_all[selected_points[0]],
     color = 'k', marker=cell_marker_i, linewidth=20, markersize=40*multiplying_factor) # start point number
    
    # plot end point
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[-1],0] + x_corr_drift_all[selected_points[-1]], rotated_cell_centers_rot[i_sel_rot,selected_points[-1],1] + y_corr_drift_all[selected_points[-1]],
     color = 'k', marker='D', linewidth=20, markersize=130) # end point black edge
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[-1],0] + x_corr_drift_all[selected_points[-1]], rotated_cell_centers_rot[i_sel_rot,selected_points[-1],1] + y_corr_drift_all[selected_points[-1]],
     color = color_i, marker='D', linewidth=20, markersize=120) # end point colored
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[-1],0] + x_corr_drift_all[selected_points[-1]], rotated_cell_centers_rot[i_sel_rot,selected_points[-1],1] + y_corr_drift_all[selected_points[-1]],
     color = 'w', marker='D', linewidth=20, markersize=80)#70) # end point white
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[-1],0] + x_corr_drift_all[selected_points[-1]], rotated_cell_centers_rot[i_sel_rot,selected_points[-1],1] + y_corr_drift_all[selected_points[-1]],
     color = 'k', marker=cell_marker_i,linewidth=20, markersize=40*multiplying_factor) # end point number

plt.imshow(imgGray, cmap='gray')
#plt.savefig("Mosaic_corrected/Figure_in_silico_mosaic_corrected_cdrift_crotation_12trajectories.png", bbox_inches='tight',pad_inches = 0)

plt.show()

1b. OR store x- and y-coordinates of all cells (200) for selected MCS (16) in array¶

In [39]:
# Define number of cells and the selected time points you want to plot the trajectory for
No_cells = 200 # number of cells that are being tracked

# Define list of arrays for x- and y-coordinates for each tracked cell
path_x_sim_all = [None] * (No_cells)
path_y_sim_all = [None] * (No_cells)
path_x_sim_all_final = [None] * No_cells
path_y_sim_all_final = [None] * No_cells
In [40]:
# Store x- and y-coordinates for each tracked cell in lists 'path_x_sim_all_final' and 'path_y_sim_all_final'
for i in range(1,int(No_cells)+1,1):
    filepath = ""
    file_dir = ['Mosaic_corrected/cell_','.txt']
    filepath = file_dir[0]+str(i)+file_dir[1]
    with open(filepath) as f:
            for line in f:
                currentline = line.split(",")
                x_line = float(currentline[0])
                y_line = float(currentline[1])
                ii = int(i)-1
                path_x_sim_all[ii] = np.append(path_x_sim_all[ii],x_line)
                path_y_sim_all[ii] = np.append(path_y_sim_all[ii],y_line)
    path_x_sim_all_final[ii] = path_x_sim_all[ii][1:(len(path_x_sim_all[ii])+1)]
    path_y_sim_all_final[ii] = path_y_sim_all[ii][1:(len(path_y_sim_all[ii])+1)]
In [41]:
# Store x- and y-coordinates of each tracked cell for only the selected time points in lists 'x_sel_final' and 'y_sel_final'
x_sel_all = [None] * No_cells
y_sel_all = [None] * No_cells
x_sel_all_final = [None] * No_cells
y_sel_all_final = [None] * No_cells

for j in range(No_cells):
    for k in range(len(selected_points)):
        x_sel_all[j] = np.append(x_sel_all[j],2*path_x_sim_all_final[j][selected_points[k]]) # this coding can probably be improved 
        y_sel_all[j] = np.append(y_sel_all[j],2*path_y_sim_all_final[j][selected_points[k]])
    x_sel_all_final[j] = x_sel_all[j][1:(len(selected_points)+1)]
    y_sel_all_final[j] = y_sel_all[j][1:(len(selected_points)+1)]

6. Plot: direction vector of each cell on final image using rotated-rotated cell centers, and correction for drift¶

In [72]:
# Color indicates orientation of direction vector
import matplotlib.colors as colors
import matplotlib.cm as cmx
import matplotlib as mpl
name = 'twilight'# or 'hsv'
cmap=mpl.colormaps[name]
#cmap = plt.cm.jet
cNorm  = colors.Normalize(vmin=-np.pi, vmax=np.pi)
scalarMap = cmx.ScalarMappable(norm=cNorm,cmap=cmap)

%matplotlib inline
imgplot_999= plt.imread('Mosaic_corrected/image0099900.png') 
plt.rcParams['figure.figsize'] = [200, 100]

# Make image gray (so no dark blue cells)
R_999, G_999, B_999 = imgplot_999[:,:,0], imgplot_999[:,:,1], imgplot_999[:,:,2]
imgGray_999 = 0.1 * R_999 + 0.1 * G_999 + 0.7*B_999

for i in range(No_cells):
    
    # Rotated and rotated back, vector downscaled by factor 20
    x_dx_rot_rot_drift = ((rotated_cell_centers_rot[i,t_end,0] - rotated_cell_centers_rot[i,t_start,0]) - (xc[t_end] - xc[t_start]))/20
    y_dy_rot_rot_drift = ((rotated_cell_centers_rot[i,t_end,1] - rotated_cell_centers_rot[i,t_start,1]) - (yc[t_end] - yc[t_start]))/20
    
    x_arrow_base_rot_rot = rotated_cell_centers_rot[i,t_end,0] - x_dx_rot_rot_drift
    y_arrow_base_rot_rot = rotated_cell_centers_rot[i,t_end,1] - y_dy_rot_rot_drift
    
    angle_cell = np.arctan2((y_dy_rot_rot_drift),(x_dx_rot_rot_drift))
    colorVal = scalarMap.to_rgba(angle_cell)
    plt.arrow(x_arrow_base_rot_rot, y_arrow_base_rot_rot, x_dx_rot_rot_drift, y_dy_rot_rot_drift, width = 1, head_width=6, head_length=6, color=colorVal)
    plt.arrow(x_arrow_base_rot_rot+300, y_arrow_base_rot_rot-290, x_dx_rot_rot_drift, y_dy_rot_rot_drift, width = 1, head_width=6, head_length=6, color=colorVal)
    
# Save figure
plt.imshow(imgGray_999, cmap='gray', alpha = .5)
#plt.savefig("Mosaic_corrected/Figure_in_silico_mosaic_corrected_t999_all_vectors_corr_drift_rot_DIV20_angle_twilight.png", bbox_inches='tight',pad_inches = 0)
plt.show()

# https://stackoverflow.com/questions/18748328/plotting-arrows-with-different-color-in-matplotlib
# use https://stackoverflow.com/questions/19576495/color-matplotlib-quiver-field-according-to-magnitude-and-direction
In [74]:
# Plot color wheel
%matplotlib inline
from matplotlib import cm
import matplotlib as mpl

fig = plt.figure()

display_axes = fig.add_axes([0.1,0.1,0.8,0.8], projection='polar')
norm = mpl.colors.Normalize(0,2*np.pi)
quant_steps = 2056
cb = mpl.colorbar.ColorbarBase(display_axes, cmap=cm.get_cmap('twilight_shifted',quant_steps),norm=norm, orientation='horizontal')

# aesthetics - get rid of border and axis labels                                   
cb.outline.set_visible(False)                                 
display_axes.set_axis_off()
plt.savefig("Mosaic_corrected/color_wheel_twilight.png", bbox_inches='tight',pad_inches = 0)
plt.show()
/var/folders/9w/72zthbm57f33y8wq5zqg5kjw0000gn/T/ipykernel_1167/1637029653.py:12: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  cb = mpl.colorbar.ColorbarBase(display_axes, cmap=cm.get_cmap('twilight_shifted',quant_steps),norm=norm, orientation='horizontal')

7a. Plot: 'Movement along long axis' vs 'End position projected on long axis', normalized¶

In [75]:
# Determine the length of the simulated gastruloid at t_end = 2 * long arm of ellipse fit. The long arm should be stored in the
# parameter 'b', and the short arm in parameter 'a', but the fitting tool sometimes stored them incorrectly.
# Therefore, check the value of a and b first to determine which parameter corresponds to the long or short arm.
# In our in silico example, the long arm of the fitted ellipses is on the x-axis. The short arm on the y-axis. 
# If a and b are incorrectly switched, the values stored in the y coordinates correspond to the long arm, and x coordinates correspond to the short arm.

if a[t_end] > b[t_end]:
    print('a > b, so x-axis is long axis')
    length_gastruloid = 2*a[t_end] # used to normalize the position of the cells along the long axis
    
    end_pos_on_long_axis = rotated_cell_centers[:,t_end,0] # rotated x-coordinate at t = t_end
    end_pos_on_long_axis_norm = (end_pos_on_long_axis - xc[t_end])/(length_gastruloid) # (rotated x-coordinate at t_end - x_center_of_ellipse at t_end) divided by length gastruloid
    
    # (rotated x-coordinate at t_end - rotated x-coordinate at t_start) - difference between center of ellipses to correct for drift
    moved_end_norm = ((end_pos_on_long_axis - rotated_cell_centers[:,t_start,0]) - (xc[t_end] - xc[t_start]) ) / length_gastruloid # divided by length of gastruloid to normalize
    
else:
    print('a < b, so y-axis is long axis')
    length_gastruloid = 2*b[t_end]
    
    # Rotated cell centers around ellipse center
    end_pos_on_long_axis = rotated_cell_centers[:,t_end,1] # rotated y-coordinate at t = t_end
    end_pos_on_long_axis_norm = (end_pos_on_long_axis - yc[t_end])/(length_gastruloid) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by length gastruloid
    
    # rotated y-coordinate at t_end - rotated y-coordinate at t_start
    moved_end_norm = ((end_pos_on_long_axis - rotated_cell_centers[:,t_start,1]) - (yc[t_end] - yc[t_start]) ) / length_gastruloid # to normalize, divide by length of gastruloid
    
# In vitro time lapse data
# See cells under "In vitro"

# -------- FIGURE ----------
# Figure parameters    
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'  
plt.xlabel('End position projected on long axis', fontsize = 16)
plt.ylabel('Movement along long axis', fontsize = 16)
plt.xlim(-0.55, 0.55)
plt.ylim(-0.51, 0.51)
plt.xticks([-0.5, -0.25, 0, 0.25, 0.5],['-L/2','-L/4','0','L/4','L/2'])
plt.yticks([-0.5, -0.25, 0, 0.25, 0.5],['-L/2','-L/4','0','L/4','L/2'])

# In vitro time lapse data
ax.scatter(end_pos_on_long_axis_iv_norm, moved_end_iv_norm, color = '#cc79a7', marker = 'o', label = 'in vitro',s = 50)

# Plot rotated cell centers around ellipse center
ax.scatter(end_pos_on_long_axis_norm, moved_end_norm, s = 5, color = 'k', label = 'in silico')
handles, labels = ax.get_legend_handles_labels()
lgd = ax.legend(handles, labels, loc='upper left', bbox_to_anchor=(1,1))
fig.set_size_inches((110/5)/2.54, (76/5)/2.54)

# Save figure
#plt.savefig("Mosaic_corrected/Figure_corr_movement_along_long_axis_legend.svg", bbox_extra_artists=([lgd]), bbox_inches='tight')

plt.show()
a > b, so x-axis is long axis

7b. Plot: 'Movement along short axis' vs 'End position projected on short axis', normalized, colored¶

In [77]:
# Plot movement along short axis

if a[t_end] > b[t_end]: # if a is long axis, rotated y-axis is short axis
    print('a > b, so y-axis is short axis')
    length_gastruloid_short = 2*b[t_end] # short axis
    
    # Rotated cell centers around ellipse center 
    end_pos_on_short_axis = rotated_cell_centers[:,t_end,1] # rotated y-coordinate at t = t_end
    end_pos_on_short_axis_norm = (end_pos_on_short_axis - yc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by (short) length gastruloid
    
    # (rotated y-coordinate at t_end - rotated y-coordinate at t_start) - difference between center of ellipses to correct for drift
    moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,1]) - (yc[t_end] - yc[t_start]) ) / length_gastruloid_short # divided by length of gastruloid to normalize
    
else:  # if b is long axis, rotated x-axis is short axis
    print('a < b, so x-axis is short axis')
    length_gastruloid_short = 2*a[t_end]
    
    # Rotated cell centers around ellipse center
    end_pos_on_short_axis = rotated_cell_centers[:,t_end,0] # rotated x-coordinate at t = t_end
    end_pos_on_short_axis_norm = (end_pos_on_short_axis - xc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by length gastruloid
    
    # rotated x-coordinate at t_end - rotated x-coordinate at t_start
    moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,0]) - (xc[t_end] - xc[t_start]) ) / length_gastruloid_short # to normalize, divide by length of gastruloid
    
# In vitro time lapse data
# See cells under "In vitro"

# -------- FIGURE ----------
# Figure parameters    
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'  
plt.xlabel('End position projected on short axis', fontsize = 16)
plt.ylabel('Movement along short axis', fontsize = 16)
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.xticks([-1, -0.5, 0, 0.5, 1],['S','-S/2','0','S/2', 'S'])
plt.yticks([-1,-0.5, 0, 0.5, 1],['-S','-S/2','0','S/2','S'])

ax.scatter(end_pos_on_short_axis_iv_norm, moved_end_iv_short_norm, color = '#cc79a7', marker = 'o', s = 50)#label = 'in vitro',s = 50)
ax.scatter(end_pos_on_short_axis_norm, moved_end_short_norm, s = 5, color = 'k')#, label = 'in silico')

fig.set_size_inches((110/5)/2.54, (76/5)/2.54)

# Save figure
#plt.savefig("Mosaic_corrected/Figure_corr_movement_along_short_axis_end_pos_short_eqaxS.svg", bbox_extra_artists=([lgd]), bbox_inches='tight')

plt.show()
a > b, so y-axis is short axis
In [78]:
# Plot movement along short axis

if a[t_end] > b[t_end]: # if a is long axis, rotated y-axis is short axis
    print('a > b, so y-axis is short axis')
    length_gastruloid_short = 2*b[t_end] # short axis
    
    # Rotated cell centers around ellipse center 
    end_pos_on_short_axis = rotated_cell_centers[:,t_end,1] # rotated y-coordinate at t = t_end
    end_pos_on_short_axis_norm = (end_pos_on_short_axis - yc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by (short) length gastruloid
    
    # (rotated y-coordinate at t_end - rotated y-coordinate at t_start) - difference between center of ellipses to correct for drift
    moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,1]) - (yc[t_end] - yc[t_start]) ) / length_gastruloid_short # divided by length of gastruloid to normalize
    
else:  # if b is long axis, rotated x-axis is short axis
    print('a < b, so x-axis is short axis')
    length_gastruloid_short = 2*a[t_end]
    
    # Rotated cell centers around ellipse center
    end_pos_on_short_axis = rotated_cell_centers[:,t_end,0] # rotated x-coordinate at t = t_end
    end_pos_on_short_axis_norm = (end_pos_on_short_axis - xc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by length gastruloid
    
    # rotated x-coordinate at t_end - rotated x-coordinate at t_start
    moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,0]) - (xc[t_end] - xc[t_start]) ) / length_gastruloid_short # to normalize, divide by length of gastruloid
    
inw_outw = np.arange(0)

for i in range(0,len(end_pos_on_short_axis_norm)):
    #print(i)
    if end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] < 0:
        inw_outw = np.append(inw_outw, 1) # inwards
    elif end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] > 0 and moved_end_short_norm[i] > end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 2) # inwards, outwards
    elif end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] > 0 and moved_end_short_norm[i] < end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 0) # outwards
    elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] < 0 and moved_end_short_norm[i] < end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 2) # inwards, outwards
    elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] < 0 and moved_end_short_norm[i] > end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 0) # outwards
    elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] > 0:
        inw_outw = np.append(inw_outw, 1) # inwards
    else:
        inw_outw = np.append(inw_outw, 3) # something is zero

# In vitro time lapse data
# See cells under "In vitro"

# -------- FIGURE ----------
# Figure parameters    
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'  
plt.xlabel('End position projected on short axis', fontsize = 16)
plt.ylabel('Movement along short axis', fontsize = 16)
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.xticks([-1, -0.5, 0, 0.5, 1],['S','-S/2','0','S/2', 'S'])
plt.yticks([-1,-0.5, 0, 0.5, 1],['-S','-S/2','0','S/2','S'])

# In vitro time lapse data
inw_outw_iv = np.arange(0)
for i in range(0,len(end_pos_on_short_axis_iv_norm)):
    #print(i)
    if end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] < 0:
        inw_outw_iv = np.append(inw_outw_iv, 1) # inwards
    elif end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] > 0 and moved_end_iv_short_norm[i] > end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 2) # inwards, outwards
    elif end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] > 0 and moved_end_iv_short_norm[i] < end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 0) # outwards
    elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] < 0 and moved_end_iv_short_norm[i] < end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 2) # inwards, outwards
    elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] < 0 and moved_end_iv_short_norm[i] > end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 0) # outwards
    elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] > 0:
        inw_outw_iv = np.append(inw_outw_iv, 1) # inwards
    else:
        inw_outw_iv = np.append(inw_outw_iv, 3) # something is zero


color_i = ['blue', #outwards
           'orange', #inwards
           'skyblue', #inwards outwards
           'yellow'] # not sure]

for i in range(0,len(end_pos_on_long_axis_iv_norm)):
    ax.scatter(end_pos_on_short_axis_iv_norm[i], moved_end_iv_short_norm[i], color = color_i[inw_outw_iv[i]], marker = 'o',s = 50)

for i in range(0,len(end_pos_on_short_axis_norm)):
    ax.scatter(end_pos_on_short_axis_norm[i], moved_end_short_norm[i], s = 5, color = color_i[inw_outw[i]])#, label = 'in silico')

fig.set_size_inches((110/5)/2.54, (76/5)/2.54)

# Save figure
#plt.savefig("Mosaic_corrected/Figure_corr_movement_along_short_axis_end_pos_short_colored_eqaxS.svg", bbox_extra_artists=([lgd]), bbox_inches='tight')

plt.show()
a > b, so y-axis is short axis

7c. Plot: 'Movement along short axis' vs. 'End position projected on long axis', normalized, colored¶

In [79]:
# Plot movement along short axis vs. end position projected on long axis

if a[t_end] > b[t_end]: # if a is long axis, rotated y-axis is short axis
    print('a > b, so y-axis is short axis')
    length_gastruloid_short = 2*b[t_end] # short axis
    length_gastruloid_long = 2*a[t_end]
    # Rotated cell centers around ellipse center 
    end_pos_on_short_axis = rotated_cell_centers[:,t_end,1] # rotated y-coordinate at t = t_end
    end_pos_on_long_axis = rotated_cell_centers[:,t_end,0] # rotated x-coordinate at t = t_end
    end_pos_on_short_axis_norm = (end_pos_on_short_axis - yc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by (short) length gastruloid
    end_pos_on_long_axis_norm = (end_pos_on_long_axis - xc[t_end])/(length_gastruloid_long) # (rotated x-coordinate at t_end - y_center_of_ellipse at t_end) divided by (short) length gastruloid
    
    # (rotated y-coordinate at t_end - rotated y-coordinate at t_start) - difference between center of ellipses to correct for drift
    moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,1]) - (yc[t_end] - yc[t_start]) ) / length_gastruloid_short # divided by length of gastruloid to normalize
    moved_end_long_norm = ((end_pos_on_long_axis - rotated_cell_centers[:,t_start,0]) - (xc[t_end] - xc[t_start]) ) / length_gastruloid_long # divided by length of gastruloid to normalize

else:  # if b is long axis, rotated x-axis is short axis
    print('a < b, so x-axis is short axis, y-axis is long axis')
    length_gastruloid_short = 2*a[t_end]
    length_gastruloid_long = 2*b[t_end]
    
    # Rotated cell centers around ellipse center
    end_pos_on_short_axis = rotated_cell_centers[:,t_end,0] # rotated x-coordinate at t = t_end
    end_pos_on_long_axis = rotated_cell_centers[:,t_end,1] # rotated y-coordinate at t = t_end
    end_pos_on_short_axis_norm = (end_pos_on_short_axis - xc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by length gastruloid
    end_pos_on_long_axis_norm = (end_pos_on_long_axis - yc[t_end])/(length_gastruloid_long) # (rotated x-coordinate at t_end - y_center_of_ellipse at t_end) divided by (short) length gastruloid
    
    # rotated x-coordinate at t_end - rotated x-coordinate at t_start
    moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,0]) - (xc[t_end] - xc[t_start]) ) / length_gastruloid_short # to normalize, divide by length of gastruloid
    moved_end_long_norm = ((end_pos_on_long_axis - rotated_cell_centers[:,t_start,1]) - (yc[t_end] - yc[t_start]) ) / length_gastruloid_long # divided by length of gastruloid to normalize

inw_outw = np.arange(0)

for i in range(0,len(end_pos_on_short_axis_norm)):
    #print(i)
    if end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] < 0:
        inw_outw = np.append(inw_outw, 1) # inwards
    elif end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] > 0 and moved_end_short_norm[i] > end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 2) # inwards, outwards
    elif end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] > 0 and moved_end_short_norm[i] < end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 0) # outwards
    elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] < 0 and moved_end_short_norm[i] < end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 2) # inwards, outwards
    elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] < 0 and moved_end_short_norm[i] > end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 0) # outwards
    elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] > 0:
        inw_outw = np.append(inw_outw, 1) # inwards
    else:
        inw_outw = np.append(inw_outw, 3) # something is zero
    
# In vitro time lapse data
# See cells under "In vitro"
inw_outw_iv = np.arange(0)
for i in range(0,len(end_pos_on_short_axis_iv_norm)):
    #print(i)
    if end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] < 0:
        inw_outw_iv = np.append(inw_outw_iv, 1) # inwards
    elif end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] > 0 and moved_end_iv_short_norm[i] > end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 2) # inwards, outwards
    elif end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] > 0 and moved_end_iv_short_norm[i] < end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 0) # outwards
    elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] < 0 and moved_end_iv_short_norm[i] < end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 2) # inwards, outwards
    elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] < 0 and moved_end_iv_short_norm[i] > end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 0) # outwards
    elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] > 0:
        inw_outw_iv = np.append(inw_outw_iv, 1) # inwards
    else:
        inw_outw_iv = np.append(inw_outw_iv, 3) # something is zero

# -------- FIGURE ----------
# Figure parameters    
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'  
plt.xlabel('End position projected on long axis', fontsize = 16)
plt.ylabel('Movement along short axis', fontsize = 16)
plt.xlim(-0.55, 0.55)
plt.ylim(-1, 1)#-0.38, 0.38)
plt.xticks([-0.5, -0.25, 0, 0.25, 0.5],['-L/2','-L/4','0','L/4','L/2'])
plt.yticks([-1,-0.5, 0, 0.5,1],['-S','-S/2','0','S/2','S'])

color_i = ['blue', #outwards
           'orange', #inwards
           'skyblue', #inwards outwards
           'yellow'] # not sure]

# In vitro time lapse data
for i in range(0,len(end_pos_on_long_axis_iv_norm)):
    ax.scatter(end_pos_on_long_axis_iv_norm[i], moved_end_iv_short_norm[i], color = color_i[inw_outw_iv[i]], marker = 'o',s = 50)

# Plot rotated cell centers around ellipse center
for i in range(0,len(end_pos_on_short_axis_norm)):
    ax.scatter(end_pos_on_long_axis_norm[i], moved_end_short_norm[i], s = 5, color = color_i[inw_outw[i]])#, label = 'in silico')

fig.set_size_inches((110/5)/2.54, (76/5)/2.54)

# Save figure
#plt.savefig("Mosaic_corrected/Figure_corr_movement_along_short_axis_end_pos_long_colored_more_spaceS.svg", bbox_extra_artists=([lgd]), bbox_inches='tight')

plt.show()
a > b, so y-axis is short axis

Noise2Void comparison¶

Calculate 5th percentile of both raw and predicted image data => use that value as the minimum displayed value¶

In [80]:
# First, change image type from 16 bit (raw tiff file) or 32 bit (prediction) to 8 bit, e.g. in FIJI.

# Read RGB images
%matplotlib inline

# Read raw and predicted image
#image_raw = cv2.imread('Mosaic_TL_mCherry/N2V_comp/T1-z11-raw-crop_8bit.tif',0)
#image_pred = cv2.imread('Mosaic_TL_mCherry/N2V_comp/prediction-t1-200epochs-z11_8bit_crop.tif',0)

#image_raw = cv2.imread('Mosaic_TL_mCherry/N2V_comp/t5-z10-crop_8bit.tif',0)
#image_pred = cv2.imread('Mosaic_TL_mCherry/N2V_comp/prediction-t5_t7-model-200epochs-z10-crop_8bit.tif',0)

#image_raw = cv2.imread('Mosaic_TL_mCherry/N2V_comp/t23-z6-crop_8bit.tif',0)
#image_pred = cv2.imread('Mosaic_TL_mCherry/N2V_comp/prediction-t23_t31-model-200epochs_z6_crop_8bit.tif',0)

image_raw = cv2.imread('Mosaic_TL_mCherry/N2V_comp/MAX_T31_crop-raw_8bit.tif',0)
image_pred = cv2.imread('Mosaic_TL_mCherry/N2V_comp/MAX_prediction-t31-200epochs-crop_8bit.tif',0)

plt.imshow(image_raw,cmap='gray')
plt.show()
plt.imshow(image_pred,cmap='gray')
plt.show()
In [81]:
# create the histogram
counts_raw, bins_raw = np.histogram(image_raw, bins=256, range=(0, 255))
counts_pred, bins_pred = np.histogram(image_pred, bins=256, range=(0, 255))

%matplotlib inline
plt.bar(bins_raw[:-1],counts_raw,width=1)
plt.bar(bins_pred[:-1],counts_pred,width=1)
Out[81]:
<BarContainer object of 256 artists>
In [82]:
# flatten the 2D image arrays
image_raw_flat = image_raw.flatten()
image_pred_flat = image_pred.flatten()

# Calculate 5th percentile of raw data, such that that value can be used as the minimum displayed value
image_raw_sort = np.sort(image_raw_flat)
print(image_raw_sort)
q5_raw=int(np.round(0.05*(len(image_raw_flat))))
print(q5_raw)
print('minimum displayed value for raw 8-bit image is: ', image_raw_sort[q5_raw])

# Calculate 5th percentile of prediction data, such that that value can be used as the minimum displayed value
image_pred_sort = np.sort(image_pred_flat)
print(image_pred_sort)
q5_pred=int(np.round(0.05*(len(image_pred_flat))))
print(q5_pred)
print('minimum displayed value for pred. 8-bit image is: ', image_pred_sort[q5_pred])
[  2   3   3 ... 255 255 255]
22714
minimum displayed value for raw 8-bit image is:  16
[  1   1   1 ... 255 255 255]
22714
minimum displayed value for pred. 8-bit image is:  4
In [ ]: